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We present a new three-dimensional general relativistic hydrodynamics code which is intended for 
simulations of stellar core collapse to a neutron star, as well as pulsations and instabilities of rotating 
relativistic stars. Contrary to the common approach followed in most existing three-dimensional 
numerical relativity codes which are based in Cartesian coordinates, in this code both the metric 
and the hydrodynamics equations are formulated and solved numerically using spherical polar coor- 
dinates. A distinctive feature of this new code is the combination of two types of accurate numerical 
schemes specifically designed to solve each system of equations. More precisely, the code uses spectral 
methods for solving the gravitational field equations, which are formulated under the assumption of 
the conformal flatness condition (CFC) for the three-metric. Correspondingly, the hydrodynamics 
equations are solved by a class of finite difference methods called high-resolution shock-capturing 
schemes, based upon state-of-the-art Riemann solvers and third-order cell-reconstruction proce- 
dures. We demonstrate that the combination of a finite difference grid and a spectral grid, on 
which the hydrodynamics and metric equations are respectively solved, can be successfully accom- 
plished. This approach, which we call Mariage des Maillages (French for grid wedding), results in 
high accuracy of the metric solver and, in practice, allows for fully three-dimensional applications 
using computationally affordable resources, along with ensuring long term numerical stability of the 
evolution. We compare our new approach to two other, finite difference based, methods to solve the 
metric equations which we already employed in earlier axisymmetric simulations of core collapse. A 
variety of tests in two and three dimensions is presented, involving highly perturbed neutron star 
spacetimes and (axisymmetric) stellar core collapse, which demonstrate the ability of the code to 
handle spacetimes with and without symmetries in strong gravity. These tests are also employed to 
assess the gravitational waveform extraction capabilities of the code, which is based on the Newto- 
nian quadrupole formula. The code presented here is not limited to approximations of the Einstein 
equations such as CFC, but it is also well suited, in principle, to recent constrained formulations of 
the metric equations where elliptic equations have a preeminence over hyperbolic equations. 

PACS numbers: 04.25.Dm, 04.30.Db, 97.60.Bw, 02.70.Bf, 02.70.Hm 



I. INTRODUCTION 

A. Relativistic core collapse simulations 

Improving our understanding of the formation of neu- 
tron stars as a result of the gravitational collapse of the 
core of massive stars is a difficult endeavour involving 
many aspects of extreme and not very well understood 
physics of the supernova explosion mechanism [lj. Nu- 
merical simulations of core collapse supernova are driv- 
ing progress in the field despite the limited knowledge 
on issues such as realistic precollapse stellar models (in- 
cluding rotation) or realistic equation of state, as well as 
numerical limitations due to Boltzmann neutrino trans- 
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port, multidimensional hydrodynamics, and relativistic 
gravity. Axisymmetric and three-dimensional approaches 
based on Newtonian gravity are available since a few 
decades now (see e.g. |2| and references therein). These 
approaches, which are constantly improving over time, 
have provided valuable information on important issues 
such as the dynamics of the collapse of a stellar core to 
nuclear density, the formation of a proto-neutron star, 
and the propagation of the shock front which ultimately 
is believed to eject the outer layers of the stellar progen- 
itor. Currently, however, even the most realistic simula- 
tions of both nonrotating and rotating progenitor models 
do not succeed in producing explosions (see Q and ref- 
erences therein). 

In addition, the incorporation of full relativistic grav- 
ity in the simulations is likely to bring in well-known 
difficulties of numerical relativity, where the attempts 
are traditionally hampered by challenging mathematical, 
computational, and algorithmic issues as diverse as the 
formulation of the field equations, robustness, efficiency, 
and long-term stability (particularly if curvature singu- 
larities are either initially present or develop during black 
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hole formation). As high densities and velocities are in- 
volved in combination with strong gravitational fields, 
gravitational collapse and neutron star formation consti- 
tute a challenging problem for general relativistic hydro- 
dynamic simulations. The pace of the progress is, no 
wonder, slow; for instance, in the three-dimensional case, 
there is still no description of core collapse in full general 
relativity today, even for the simplest matter models one 
can conceive, where all microphysics is neglected. 

In recent years, the interest in performing core collapse 
simulations has been further motivated by the necessity 
of obtaining reliable gravitational waveforms from (ro- 
tating) core collapse, one of the main targets of gravi- 
tational radiation for the present and planned interfer- 
ometer detectors such as LIGO, GEO600, and VIRGO 
(see for a review). As a result of the complexities 
listed above, it is not surprising that most previous stud- 
ies aimed at computing the gravitational wave signature 
of core collapse supernovae have considered greatly sim- 
plified parameterized models 0, H B II II U E E E 
E Q Hi E EE In addition to the burst signal 
of gravitational waves emitted during core bounce, mul- 
tidimensional simulations have also provided the signals 
produced by convection El ( see a l so E f° r f ne m °st re- 
alistic simulations available at present), as well as those 
from the resulting neutrino emission Ei |^| . 

From the above references it becomes apparent that 
our understanding of core collapse and neutron star 
formation has advanced mainly by studies carried out 
employing Newtonian dynamics. The situation is now 
slowly changing, at least for simplified matter models 
where microphysics and radiation transport are not yet 
included, with new formulations of the Einstein field 
equations and of the general relativistic hydrodynamics 
equations. Unfortunately, the 3+1 Einstein equations de- 
scribing the dynamics of spacetime are a complicated set 
of coupled, highly nonlinear hyperbolic-elliptic equations 
with plenty of terms. Their formulation in a form suit- 
able for accurate and stable numerical calculations is not 
unique, and constitutes one of the major fields of current 
research in numerical relativity (see E, and refer- 
ences therein) . Not surprisingly, approximations of those 
equations have been suggested, such as the conformal 
flatness condition of Isenberg- Wilson-Mathews E, E 
(CFC hereafter), who proposed to approximate the 3- 
metric of the 3 + 1 decomposition by a conformally flat 
metric. 

Using this approximation, Dimmelmeier et al. [Tol ITU 
El presented the first relativistic simulations of the core 
collapse of rotating polytropes and neutron star forma- 
tion in axisymmetry, providing an in-depth analysis of 
the dynamics of the process as well as of the gravita- 
tional wave emission. The results showed that relativis- 
tic effects may qualitatively change in some cases the 
dynamics of the collapse obtained in previous Newto- 
nian simulations In particular, core collapse with 
multiple bounces was found to be strongly suppressed 
when employing relativistic gravity. In most cases, com- 



pared to Newtonian simulations, the gravitational wave 
signals are weaker and their spectra exhibit higher av- 
erage frequencies, as the newly born proto-neutron stars 
have stronger compactness in the deeper relativistic grav- 
itational potential. Therefore, telling from simulations 
based on rotating polytropes, the prospects for detec- 
tion of gravitational wave signals from supernovae are 
most likely not enhanced by taking into account relativis- 
tic gravity. The gravitational wave signals computed by 
Dimmelmeier et al. E> E> 03 are within the sensitivity 
range of the planned laser interferometer detectors if the 
source is located within our Galaxy or in its local neigh- 
bourhood. A catalogue of the core collapse waveforms 
presented in El is available electronically E ■ This cat- 
alogue is currently being employed by gravitational wave 
data analysis groups to calibrate their search algorithms 
(see e.g. [27j for results concerning the VIRGO group). 

More recently, Shibata and Sekiguchi [17J have pre- 
sented simulations of axisymmetric core collapse of rotat- 
ing polytropes to neutron stars in full general relativity. 
These authors used a conformal-traceless reformulation 
of the 3+1 gravitational field equations commonly re- 
ferred to in the literature by the acronym BSSN after 
the works of E> E (but note that many of the new fea- 
tures of the BSSN formulation were anticipated as early 
as 1987 by Nakamura, Oohara, and Kojima |3(j)- The 
results obtained for initial models similar to those of El 
agree to high precision in both the dynamics of the col- 
lapse and the gravitational waveforms. This conclusion, 
in turn, implies that, at least for core collapse simulations 
to neutron stars, CFC is a very precise approximation of 
general relativity. 

We note that in the relativistic core collapse simula- 
tions mentioned thus far Ej Ej the gravitational ra- 
diation is computed using the (Newtonian) quadrupole 
formalism. To the best of our knowledge the only ex- 
ception to this is the work of Siebel et al. [3l]|, where, 
owing to the use of the characteristic (light-cone) formu- 
lation of the Einstein equations, the gravitational radi- 
ation from axisymmetric core collapse simulations was 
unambiguously extracted at future null infinity without 
any approximation. 



B. Einstein equations and spectral methods 

The most common approach to numerically solve 
the Einstein equations is by means of finite differences 
(see [23 and references therein). However, it is well 
known that spectral methods E, El are ^ ar more accu- 
rate than finite differences for smooth solutions (e.g. best 
for initial data without discontinuities), being particu- 
larly well suited to solve elliptic and parabolic equations. 
Good results can be obtained for hyperbolic equations 
as well, as long as no discontinuities appear in the solu- 
tion. The basic principle underlying spectral methods is 
the representation of a given function f(x) by its coeffi- 
cients in a complete basis of orthonormal functions: sines 
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and cosines (Fourier expansion) or a family of orthogonal 
polynomials (e.g. Chebyshev polynomials T^x) or Leg- 
endre polynomials). In practice, of course, only a finite 
set of coefficients is used and one approximates / by the 
truncated series f(x) ~ J2i=o c i^i{ x ) °f such functions. 
The use of spectral methods results in a very high accu- 
racy, since the error made by this truncation decreases 
like e~ n for smooth functions (exponential convergence). 

In an astrophysical context spectral methods have al- 
lowed to study subtle phenomena such as the develop- 
ment of physical instabilities leading to gravitational col- 
lapse 34]. In the last few years, spectral methods have 
been successfully employed by the Meudon group [sfij in a 
number of relativistic astrophysics scenarios |3J| , among 
them the gravitational collapse of a neutron star to a 
black hole, the infall phase of a tri-axial stellar core in 
a core collapse supernova (extracting the gravitational 
waves emitted in such process), the construction of equi- 
librium configurations of rapidly rotating neutron stars 
endowed with magnetic fields, or the tidal interaction 
of a star with a massive black hole. Their most recent 
work concerns the computation of the inertial modes of 
rotating stars |37l ] , of quasi-equilibrium configurations of 
co-rotating binary black holes in general relativity [38| . 
as well as the evolution of pure gravitational wave space- 
times 39] . To carry out these numerical simulations the 
group has developed a fully object-oriented library called 
Lorene j4jj (based on the C++ computer language) 
to implement spectral methods in spherical coordinates. 
Spectral methods are now employed in numerical relativ- 
ity by other groups as well |4l|, |43 • 



C. Hydrodynamics equations and HRSC schemes 

On the other hand, robust finite difference schemes to 
solve hyperbolic systems of conservation (and balance) 
laws, such as the Euler equations of fluid dynamics, are 
known for a long time and have been employed success- 
fully in computational fluid dynamics (see e.g. 0] an d 
references therein). In particular, the so-called upwind 
high-resolution shock-capturing schemes (HRSC schemes 
hereafter) have shown their advantages over other type 
of methods even when dealing with relativistic flows with 
highly ultrarelativistic fluid speeds (see e.g. 0, Hg] and 
references therein). HRSC schemes are based on the 
mathematical information contained in the characteristic 
speeds and fields (eigenvalues and eigenvectors) of the Ja- 
cobian matrices of the system of partial differential equa- 
tions. This information is used in a fundamental way to 
build up either exact or approximate Riemann solvers to 
propagate forward in time the collection of local Riemann 
problems contained in the initial data, once these data 
are discretized on a numerical grid. These schemes have 
a number of interesting properties: (1) The convergence 
to the physical solution (i.e. the unique weak solution 
satisfying the so-called entropy condition) is guaranteed 
by simply writing the scheme in conservation form, (2) 



the discontinuities in the solution are sharply and stably 
resolved, and (3) these methods attain a high order of 
accuracy in smooth parts of the solution. 



D. Mariage des Maillages 

From the above considerations, it seems a promising 
strategy, in the case of relativistic problems where cou- 
pled systems of elliptic (for the spacetime) and hyperbolic 
(for the hydrodynamics) equations must be solved, to 
use spectral methods for the former and HRSC schemes 
for the latter (where discontinuous solutions may arise). 
Showing the feasibility of such an approach is, in fact, 
the main motivation and aim of this paper. Therefore, 
we present and assess here the capabilities of a new, 
fully three-dimensional code whose distinctive features 
are that it combines both types of numerical schemes 
and implements the field equations and the hydrody- 
namic equations using spherical coordinates. It should 
be emphasized that our Mariage des Maillages approach 
is hence best suited for formulations of the Einstein equa- 
tions which favor the appearance of elliptic equations 
against hyperboli c eq uations, i.e. either approximations 
such as CFC 0, 125) (the formulation we adopt in the 
simulations reported in this paper), higher-order post- 
Newtonian extensions |4a. or exact formulations as re- 
cently proposed by [39, 47]. The hybrid approach put 
forward here has a successful precedent in the literature; 
using such combined methods, first results were obtained 
in one-dimensional core collapse in the framework of a 
tensor-scalar theory of gravitation . 

We note that one of the main limitations of the pre- 
vious axisymmetric core collapse simulations presented 
m was the CPU time spent when solving 

the elliptic equations describing the gravitational field 
in CFC. The restriction was severe enough to prevent 
the practical extension of the investigation to the three- 
dimensional case. In that sense, spectral methods are 
again particularly appropriate as they provide accurate 
results with reasonable sampling, as compared with finite 
difference methods. 

The three-dimensional code we present in this paper 
has been designed with the aim of studying general rel- 
ativistic astrophysical scenarios such as rotational core 
collapse to neutron stars (and, eventually, to black holes), 
as well as pulsations and instabilities of the formed com- 
pact objects. Core collapse may involve, obviously, mat- 
ter fields which are not rotationally symmetric. While 
during the infall phase of the collapse the deviations from 
axisymmetry should be rather small, for rapidly rotating 
neutron stars which form as a result of the collapse, or 
which may be spun up by accretion at later times, ro- 
tational (nonaxisymmetric) bar mode instabilities may 
develop, particularly in relativistic gravity and for differ- 
ential rotation. In this regard, in the previous axisym- 
metric simulations of Dimmelmeier et al. |l2j , some of the 
most extremely rotating initial models yielded compact 
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remnants which are above the thresholds for the develop- 
ment of such bar mode instabilities on secular or even dy- 
namic time scales for Maclaurin spheroids in Newtonian 
gravity (which are (3 S ~ 0.14 and (3d ~ 0.27, respectively, 
with (3 = -E r /|i?b| being the ratio of rotational energy 
and gravitational binding energy). 

Presently, only a few groups worldwide have developed 
finite difference, three-dimensional (Cartesian) codes ca- 
pable of performing the kind of simulations we aim at, 
where the joint integration of the Einstein and hydro- 
dynamics equations is required 0, 15(3, • Further 3D 
codes are currently being developed by a group in the 
U.S. an( i by a E.U. Research Training Network col- 
laboration 153.1541. 



E. Organization of the paper 

The paper is organized as follows: In Section [n] we in- 
troduce the assumptions of the adopted physical model 
and the equations governing the dynamics of a general 
relativistic fluid and the gravitational held. Section IIIII 
is devoted to describing algorithmic and numerical fea- 
tures of the code, such as the setup of both the spectral 
and the finite difference grids, as well as the basic ideas 
behind the HRSC schemes we have implemented to solve 
the hydrodynamics equations. In addition, a detailed 
comparison of the three different solvers for the metric 
equations and their practical applicability is given. In 
Section IIVI we present a variety of tests of the numeri- 
cal code, comparing the metric solver based on spectral 
methods to two other alternative methods using hnite dif- 
ferences. We conclude the paper with a summary and an 
outlook to future applications of the code in Section [3 
We use a spacelike signature (—,+,+,+) and units in 
which c = G = 1 (unless explicitly stated otherwise). 
Greek indices run from to 3, Latin indices from 1 to 3, 
and we adopt the standard convention for the summation 
over repeated indices. 



specific enthalpy, defined as ft, = 1 + e + P/p, where e 
is the specific internal energy. The three-velocity of the 
fluid, as measured by an Eulerian observer at rest in a 
spacelike hypersurface Et is given by 



a 



(2) 



where a is the lapse function and f3 l is the shift vector 
(see Section lrTB|) . 

Following the work laid out in |55j we now introduce 
the following set of conserved variables in terms of the 
primitive (physical) hydrodynamic variables (p, Vi,e): 

D = pW, 
Si = phW 2 Vi , 
t = phW 2 -P - D. 

In the above expressions W is the Lorentz factor de- 
fined as W = au°, which satisfies the relation W = 
l/Vl — ViV % and w = jijV^ , where jij is the 3-metric. 

Using the above variables, the local conservation 
laws Q can be written as a first-order, flux-conservative 
hyperbolic system of equations, 



dt 



dx l 



Q, 



(3) 



with the state vector, flux vector, and source vector given 
by 

U=[D,S j ,r], 

F l = [Dv\ Sjff + 6}P, t& + Pv'] , 

ufdguj n A A ( 4 ) 



Q 



d In a 



dx^ 



1 y. v 9\j 



II. PHYSICAL MODEL AND EQUATIONS 
A. General relativistic hydrodynamics 

1. Flux-conservative hyperbolic formulation 

Let p denote the rest-mass density of the fluid, 
its four-velocity, and P its pressure. The hydrody- 
namic evolution of a relativistic perfect fluid with rest- 
mass current J M = pu^ and energy-momentum tensor 
T^ v = phu^u v -(- Pg^ v in a (dynamic) spacetime g^ v is 
determined by a system of local conservation equations, 
which read 



: 0, 



0, 



(1) 



where denotes the covariant derivative. The quan- 
tity h appearing in the energy-momentum tensor is the 



Here v l = v l — f3 l /a, and \J—g = ct^Ff, with g = det^^) 
and 7 = det(7y) being the determinant of the 4-metric 
and 3-metric, respectively (see Section Ulrj In addi- 
tion, are the Christoffel symbols associated with g^. 



2. Equation of state 

The system of hydrodynamic equations © is closed 
by an equation of state (EoS) which relates the pressure 
to some thermodynamically independent quantities, e.g. 
P = P(p,e). As in [Til lL2t l3ll] we have implemented in 
the code a hybrid ideal gas EoS (5(|, which consists of a 
polytropic pressure contribution and a thermal pressure 
contribution, P = P p + P t h- This EoS, which despite its 
simplicity is particularly suitable for stellar core collapse 
simulations, is intended to model the degeneracy pres- 
sure of the electrons and (at supranuclear densities) the 
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pressure due to nuclear forces in the polytropic part, and 
the heating of the matter by shock waves in the thermal 
part. The hybrid EoS is constructed as follows. 

For a rotating stellar core before collapse the polytropic 
relation between the pressure and the rest mass density, 



P P = Kp\ 



(5) 



with 7 = 7 ini = 4/3 and K = 4.897 x 10 14 (in cgs units) 
is a fair approximation of the density and pressure strat- 
ification Q. 

In order to start the gravitational collapse of a con- 
figuration initially in equilibrium, the effective adiabatic 
index 7 is reduced from 7i n ; to 71 on the initial time slice. 
During the infall phase of core collapse the matter is as- 
sumed to obey a polytropic EoS JSJ, which is consistent 
with the ideal gas EoS for a compressible inviscid fluid, 
P = (7 - l)pe. 

To approximate the stiffening of the EoS for densities 
larger than nuclear matter density /9 nuc , we assume that 
the adiabatic index 7 jumps from 71 to 72 at p = p nuc . 
At core bounce a shock forms and propagates out, and 
the matter accreted through the shock is heated, i.e. its 
kinetic energy is dissipated into internal energy. This is 
reflected by a nonzero Pth = peth(7th — 1), where e t h = 
e — e p with e p = P p /[p(7 — 1)], in the post-shock region. 
We choose 7 th = 1.5. This choice describes a mixture 
of relativistic (7 — 4/3) and nonrelativistic (7 = 5/3) 
components of an ideal fluid. 

Requiring that P and e are continuous at the transition 
density p rm c, one can construct an EoS for which both the 
total pressure P and the individual contributions P p and 
Pth are continuous at p n uc, and which holds during all 
stages of the collapse: 



P = ^—^Kp^-ipi 



(7th-l)(7~7i) ^ 71 -i 
(7i-l) (72-1) Pnuc ' 



7-1 

+ (7th - l)pc. (6) 

For more details about this EoS, we refer to [Tlll5^ |. 

Our implementation of the hybrid EoS allows us to 
suppress the contribution of the thermal pressure Pth- 
In this case the EoS © analytically reduces to the poly- 
tropic relation We use this EoS, with different values 
for 7 and K, in the simulations of polytropic neutron star 
models presented below. 



to a hypersurface, (3 l is the spacelike shift three- vector 
which describes the motion of coordinates within a sur- 
face, and 7y is the spatial three-metric. 

In the 3 + 1 formalism, the Einstein equations are split 
into evolution equations for the three-metric 7y and the 
extrinsic curvature Kij, and constraint equations (the 
Hamiltonian and momentum constraints) which must be 
fulfilled at every spacelike hypersurface: 

<9 t7 y = -2aKij + Vi/3j + V^-ft, 
d t K i:j = -ViVja + a(R l3 + KKi-j - 2K lk Kf) 



(8) 



-SnafSij-^-iS^-pu)), 
= R + K 2 - A , , A ' •' - 167T/9H, 

= Vi(i<r ij - y j k) - 8nS j . 



In these equations Vj is the covariant derivative with 
respect to the three-metric 7y , Py is the corresponding 
Ricci tensor, R is the scalar curvature, and K is the trace 
of the extrinsic curvature . The matter fields appear- 
ing in the above equations, Sij, , and pn = phW 2 — P, 
are the spatial components of the stress-energy tensor, 
the three momenta, and the total energy, respectively. 

The ADM equations have been repeatedly shown over 
the years to be intrinsically numerically unstable. Re- 
cently, there have been numerous attempts to reformu- 
late above equations into forms better suited for nu- 
merical investigations (see [U H3, H3, and refer- 
ences therein). These approaches to delay or entirely 
suppress the excitation of constraint violating unstable 
modes include the BSSN reformulation of the ADM sys- 
tem |23, 12^. poll (see Section IIB|I , hyperbolic reformula- 
tions (see |58j and references therein) . or a new form with 
maximally constrained evolution [39|. In our opinion a 
consensus seems to be emerging currently in numerical 
relativity, which in general establishes that the more con- 
straints are used in the formulation of the equations the 
more numerically stable the evolution is. 



2. Conformal flatness approximation for the spatial metric 



B. Metric equations 

1. ADM metric equations 

We adopt the ADM 3 + 1 formalism |57j to foliate the 
spacetime into a set of non-intersecting spacelike hyper- 
surfaces. The line element reads 

ds 2 = -a 2 dt 2 + (dx l + f3 l dt)(dx 3 + fidt), (7) 

where a is the lapse function which describes the rate of 
advance of time along a timelike unit vector normal 



Based on the ideas of Isenberg |24| and Wilson et 
al. |25|. and as it was done in the work of Dimmelmeier 
et al. [12| , we approximate the general metric g^ v by re- 
placing the spatial three-metric 7y with the conformally 
flat three-metric, 7y = <fi 4 jij, where 7^ is the flat met- 
ric (ffij = Sij in Cartesian coordinates). In general, the 
conformal factor (f> depends on the time and space co- 
ordinates. Therefore, at all times during a numerical 
simulation we assume that all off-diagonal components 
of the three-metric are zero, and the diagonal elements 
have the common factor </> 4 . 

In CFC the following relation between the time deriva- 
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tive of the conformal factor and the shift vector holds: 



d t 4> = |v fe/ 3 fe . 

6 



(9) 



With this the expression for the extrinsic curvature be- 
comes time-independent and reads 



± (Vifr + V^A - \iijV k A . (10) 



If we employ the maximal slicing condition, K = 0, then 
in the CFC approximation the ADM equations JSJ reduce 
to a set of five coupled elliptic (Poisson-like) nonlinear 
equations for the metric components, 

f 7K 
A(o0) = 2™0 5 ph(3W 2 - 2) + 5P+-^ , (11) 

\ l07T / 

A/3 1 = Wira^S* + 2c/> w K i ^ j (jj^j - iv'Vfc/S*, 

where Vj and A are the flat space Nabla and Laplace 
operators, respectively. We note that the way of writ- 
ing the metric equations with a Laplace operator on the 
left hand side can be exploited by numerical methods 
specifically designed to solve such kind of equations (see 
Sections IIII D 21 and IIIID3I below) . 

These elliptic metric equations couple to each other 
via their right hand sides, and in case of the three equa- 
tions for the components of (3 l also via the operator A 
acting on the vector /3\ They do not contain explicit 
time derivatives, and thus the metric is calculated by 
a fully constrained approach, at the cost of neglecting 
some evolutionary degrees of freedom in the spacetime 
metric. In the astrophysical situations we plan to ex- 
plore (e.g. evolution of neutron stars or core collapse of 
massive stars), the equations are entirely dominated by 
the source terms involving the hydrodynamic quantities 
p, P, and v l , whereas the nonlinear coupling through the 
remaining, purely metric, source terms becomes only im- 
portant for strong gravity. On each time slice the metric 
is hence solely determined by the instantaneous hydro- 
dynamic state, i.e. the distribution of matter in space. 

Recently, Cerda-Duran et al. 0] have extended the 
above CFC system of equations (and the corresponding- 
core collapse simulations in CFC reported in 0) by the 
incorporation of additional degrees of freedom in the ap- 
proximation, which render the spacetime metric exact up 
to the second post-Newtonian order. Despite the exten- 
sion of the five original elliptic CFC metric equations for 
the lapse, the shift vector, and the conformal factor by 
additional equations, the final system of equations in the 
new formulation is still elliptic. Hence, the same code and 
numerical schemes employed in |l2j and in the present 
work can be used. The results obtained by Cerda-Duran 
et al. 0| for a representative subset of the core collapse 
models in show only minute differences with respect 



to the CFC results, regarding both the collapse dynam- 
ics and the gravitational waveforms. We point out that 
Shibata and Sekiguchi have recently considered ax- 
isymmctric core collapse of rotating polytropes to neu- 
tron stars in full general relativity (i.e. no approxima- 
tions) using the 3+1 BSSN formulation of the Einstein 
equations. Interestingly, the results obtained for initial 
models similar to those of |l2j agree to high precision 
in the dynamics of the collapse and on the gravitational 
waveforms, which supports the suitability and accuracy 
of the CFC approximation for simulations of relativistic 
core collapse to neutron stars (see also Section HVB4|1 . 

In addition, there has been a direct comparison be- 
tween the CFC approximation and perturbative analyti- 
cal approaches (post-Newtonian and effective-one-body) , 
which shows a very good agreement in the determination 
of the innermost stable circular orbit of a system of two 
black holes |59| . 



3. Metric equation terms with noncompact support 

In general, the right hand sides of the metric equa- 
tions l|llfl contain nonlinear source terms of noncompact 
support. For a system with an isolated matter distribu- 
tion bounded by some stellar radius r s , the source term 
of each of the metric equations for a metric quantity u 
can be split into a "hydrodynamic" term with compact 
support Sh and a purely "metric" term with noncompact 
support S m . Where no matter is present, only the metric 
term remains: 



/ S h (u) + S m (u) 
S m {u) 



for r < r s 
for r > r s 



(12) 



The source term vanishes only for K%j = and thus 
j3 l = 0, i.e. if the three-velocity vanishes and the mat- 
ter is static. As a consequence of this, only a spherically 
symmetric static matter distribution will yield a time- 
independent solution to Eq. I|12f) . which is equivalent to 
the spherically symmetric Tolman-Oppenheimer-Volkoff 
(TOV) solution of hydrostatic equilibrium. In this case 
the vacuum metric is given by the solution of a homoge- 
neous Poisson equation, u = k± + A^A", the constants k\ 
and &2 being determined by boundary values e.g. at r s . 

A time-dependent spherically symmetric matter inte- 
rior suffices to yield a nonstatic vacuum metric (u = u(t) 
everywhere). However, this is not a contradiction to 
Birkhoff's theorem, as it is purely a gauge effect. A 
transformation of the vacuum part of the metric from 
an isotropic to a Schwarzschild-like radial coordinate 
leads to the static (and not conformally flat) standard 
Schwarzschild vacuum spacetime. 

Thus, in general, the vacuum metric solution to 
Eqs. Ijlll) cannot be obtained analytically, and therefore 
(except for TOV stars) no exact boundary values can be 
imposed for <f>, a, and f3 l at some finite radius r. We note 
that this property of the metric equations is no conse- 
quence of the approximative character of conformal flat- 
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ness, as in spherical symmetry the CFC renders the exact 
ADM equations JSJ, but rather results from the choice of 
the (isotropic) radial coordinate. 



of axisymmetric rotational core collapse of the order of 
10~ 4 . This small violation can be entirely attributed to 
the interaction of the stellar matter with the artificial 
atmosphere (see Appendix IA 2} . 



III. NUMERICAL METHODS 
A. Finite difference grid 

The expressions for the hydrodynamic and metric 
quantities outlined in Section [H] are in covariant form. 
For a numerical implementation of these equations, how- 
ever, we have to choose a suitable coordinate system 
adapted to the geometry of the astrophysical situations 
intended to be simulated with the code. 

As we plan to investigate isolated systems with matter 
configurations not too strongly departing from spherical 
symmetry with a spacetime obeying asymptotic flatness, 
the formulation of the hydrodynamic and metric equa- 
tions, Eqs. (PJ) and and their numerical implementa- 
tion are based on spherical polar coordinates (t,r,9,(p). 
This coordinate choice facilitates the use of fixed grid 
refinement in form of nonequidistant radial grid spac- 
ing. Additionally, in spherical coordinates the boundary 
conditions for the system of partial differential metric 
equations are simpler to impose (at finite or infinite 
distance) on a spherical surface than on a cubic surface 
if Cartesian coordinates were used. We have found no 
evidence of numerical instabilities arising at the coordi- 
nate singularities at the origin (r = 0) or at the axis 
(9 = 0, 7r) in all simulations performed thus far with the 
code (see [6(3, EOJ for related discussions on instabilities 
in codes based upon spherical coordinates). 

Both the discretized hydrodynamic and metric quan- 
tities are located on the Eulerian finite difference grid 
at cell centers (ri,9j,tpk), where i,j, k run from 1 to 
n r ,no,n v , respectively. The angular grid zones in the 
0- and 95-direction are each equally spaced, while the ra- 
dial grid, which extends out to a finite radius r^ larger 
than the stellar radius r s , can be chosen to be equally 
or logarithmically spaced. Each cell is bounded by two 
interfaces in each coordinate direction. Values on ghost 
zone cell centers, needed to impose boundary conditions, 
are obtained with the symmetry conditions described 
in [Tl| . We further assume equatorial plane symmetry 
in all simulations presented below (the code, however, is 
not restricted to this symmetry condition). Expressions 
containing finite differences in space on this grid are cal- 
culated with second order accuracy. 

Note that the space between the surface of the star, 
the radius of which in general is angular dependent, and 
the outer boundary of the finite difference grid is filled 
with an artificial atmosphere (as done in codes similar to 
ours, see IH, |5^| ) This atmosphere obeys the poly- 
tropic EoS 0, and has a very low density such that its 
presence does not affect the dynamics of the star . 
As an example, we observe a slight violation of conserva- 
tion of rest mass and angular momentum in simulations 



B. Spectral methods and grid 

1. Spectral methods 

Our most general metric solver is based on spectral 
methods (see Section 1111 D 31) . The basic principle of 
these methods has been given in Section ll Bl Let us now 
describe some details of our implementation in the case 
of 3D functions in spherical coordinates. The interested 
reader can refer to [36| for details. A function / can be 
decomposed as follows (£ is linked with the radial coor- 
dinate r, as given below): 



fit 0,*>) = ££E ciMQYfie, ip), 



(13) 



fc=0 j=0 i=0 



where Yj(6,cp) are spherical harmonics. The angular 
part of the function can also be decomposed into a 
Fourier series, to compute angular derivatives more eas- 
ily. If / is represented by its coefficients Cijk, it is easy 
to obtain the coefficients of e.g. df /dr, Af (or the result 
of any linear differential operator applied to /) thanks 
to the properties of Chebyshev polynomials or spherical 
harmonics. For instance, to compute the coefficients of 
the radial derivative of /, we make use of the following 
recursion formula on Chebyshev polynomials: 

dT n+1 (x) n+ldT n _i(x) 

= 2(n + l)T n {x) H r : Vn > 1. 



dx 



1 



dx 



(14) 

A grid is still needed for two reasons: firstly, to calcu- 
late these coefficients through the computation of inte- 
grals, and secondly to evaluate non-linear operators (e.g. 
V/ x V/), using the values of the functions at grid points 
(in physical space). The spectral grid points, called collo- 
cation points are situated at (fj, 6j, <pk), where k run 
from 1 to n r ,ng 1 n ipi respectively. They are the nodes 
of a Gauss-Lobato quadrature used to compute the in- 
tegrals giving the spectral coefficients. The use of Fast 
Fourier Transforms (FFT) for the angular part requires 
equally spaced points in the angular directions, whereas 
a fast Chebyshev transform (also making use of FFT) 
requires that the radial grid points correspond, in £, to 
the zeros of Th r . Note that in our simulations each of the 
domains contains the same number of radial and angular 
collocation points. 

In order to be able to cover the entire space (r E 
[0, +oo]) and to handle coordinate singularities at the 
origin (r — -- 0), we use several grid domains: 

• a nucleus spanning from r = to r<i, where we set 
r = a£, with £ S [0, 1] and a being a constant (we 
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use either only even Chebyshev polynomials 22i(£), 
or only odd polynomials T2i+i(£)); 

• an arbitrary number (including zero) of shells 
bounded by the inner radius r d i and outer radius 
r cH+i! where we set r = + Pi with £ e [—1, 1] 
and «i and /3j being constants depending on the 
shell number i; 

• a compactified external domain extending from the 
outer boundary of the finite difference grid at rfd 
to radial infinity, where we set r = l/[a c (£ + 1)], 
with £ G [—1, f] and a c being a constant. 

Furthermore, we assume that the ratio f d between the 
outer boundary radii of two consecutive domains is con- 
stant, which yields the relation 



/d=( ^'' 



(15) 



where rid is the number of domains (including the nucleus 
and the external compactified domain). Thus a particu- 
lar choice of na and fixing the radius of the nucleus rd 
completely specifies the setup of the spectral grid: 



r d i = r d , 



rdi = fd X r d i-x, 



(16) 



r<in d -i 

r d rid 



Tfd, 
CO. 



The setup of the spectral grid and the associated finite 
difference grid for a typical stellar core collapse model is 
exemplified in Fig. ^ for n r — 33 grid points per spectral 
radial domain and n r — 200 finite difference grid points. 
Particularly in the central parts of the star (upper panel) 
the logarithmic radial spacing of the finite difference grid 
is obvious. While the finite difference grid ends at the 
finite radius rfd (with the exception of four ghost zones, 
which are needed for the hydrodynamic reconstruction 
scheme; see Section lill CI) , the radially compactified out- 
ermost 6th domain of the spectral grid covers the entire 
space to radial infinity (lower panel). The finite differ- 
ence grid is fixed in time, while the boundaries rdi of 
the spectral radial domains (and thus the radial colloca- 
tion points) change adaptively during the evolution (for 
details, we refer to Section flV B 3D . Note that the ra- 
dial collocation points of the spectral grid, which corre- 
spond to the roots of the Chebyshev polynomials (for the 
Gauss-Lobato quadrature) , are concentrated towards the 
domain boundaries. 

Generally speaking, in order to achieve a compara- 
ble accuracy in the representation of functions and their 
derivatives, the finite difference grid needs much more 
points than the spectral one. For example, when consid- 
ering the representation of some function like exp(— x 2 ) 



on the interval [0, 
polynomials need 



1], spectral methods using Chebyshev 
~ 30 coefficients (and grid points) to 
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FIG. 1: Radial setup of the initial spectral grid (collocation 
points are marked by plus symbols) and the time-independent 
finite difference grid (cell centers are marked by filled circles, 
separated by cell interfaces symbolized by vertical dashes) for 
a typical core collapse simulation. The upper panel shows 
the innermost 500 km containing the nucleus (ending at rd ~ 
200 km), the first shell, and a part of the second shell of the 
spectral grid. In the lower panel a part of the last regular 
shell (which is confined by the outer boundary of the finite 
difference grid at rfd ~ 2200 km) and the beginning of the 
compactified domain of the spectral grid are plotted. The 
domain boundaries are indicated by vertical dotted lines. 



reach machine double precision (10 -16 ) for the represen- 
tation of the function and 10~ 13 for the representation of 
its first derivative. For comparison, a third order scheme 
based on finite differences needs ~ 10 5 points to achieve 
the same accuracy. 



2. Communication between grids 

Passing information from the spectral grid to the fi- 
nite difference grid is technically very easy. Knowing 
the spectral coefficients of a function, this step simply 
requires the evaluation of the sum (|13fl at the finite dif- 
ference grid points. The drawback of this method, as it 
will be discussed in Section IIV Al is the computational 
time spent. In 3D this time can even be larger than the 
time spent by the spectral elliptic solver. Going from 
the finite difference grid to the spectral grid requires an 
actual interpolation, taking special care to avoid Gibbs 
phenomena that can appear in the spectral representa- 
tion of discontinuous functions. The matter terms en- 
tering in the sources of the gravitational field equations 
can be discontinuous when a shock forms. Thus, it is 
necessary to smooth or filter out high frequencies that 
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would otherwise spoil the spectral representation. This 
introduces a numerical error in the fields that should re- 
main within the overall error of the code. The important 
point to notice is that an accurate description needs not 
be achieved in the spectral representation of the sources 
(the hydrodynamic quantities are well described on the fi- 
nite difference grid) , but in that of the gravitational field, 
which is always continuous, as well as its first derivatives. 

Technically, we interpolate from the finite difference 
grid to the spectral grid using a one-dimensional algo- 
rithm and intermediate grids. We first perform an inter- 
polation in the r-direction, then in the (^-direction and 
finally in the (^-direction. We can choose between piece- 
wise linear or parabolic interpolations, and a scheme that 
globally minimizes the norm of the second derivative of 
the interpolated function j^. The filtering of spectral 
coefficients is performed a posteriori by removing the co- 
efficients corresponding to higher frequencies. For exam- 
ple, in the radial direction, this is done by canceling the 
Cijk in Eq. fTSjl for i larger than a given threshold. In 
practice, best results were found when cancelling the last 
third of radial coefficients. This can be linked with the so- 
called "two-thirds rule" used for spectral computations 
of quadratically nonlinear equations Nevertheless, 
a different (higher) threshold would also give good re- 
sults, in the sense that there are no high-frequency terms 
rising during the metric iteration. 



C. High-resolution shock-capturing schemes 

As in our previous axisymmetric code 0, ^| , in the 
present code the numerical integration of the system of 
hydrodynamic equations is performed using a Godunov- 
type scheme. Such schemes are specifically designed to 
solve nonlinear hyperbolic systems of conservation laws 
(see, e.g. f° r general definitions and for spe- 

cific details regarding their use in special and general rcla- 
tivistic hydrodynamics). In a Godunov-type method the 
knowledge of the characteristic structure of the equations 
is crucial to design a solution procedure based upon either 
exact or approximate Riemann solvers. These solvers, 
which compute at every cell-interface of the numerical 
grid the solution of local Riemann problems, guarantee 
the proper capturing of all discontinuities which may ap- 
pear in the flow. 

The time update of the hydrodynamic equations © 
from i" to t n+1 is performed using a method of lines in 
combination with a second-order (in time) conservative 
Runge-Kutta scheme. The basic conservative algorithm 
reads: 

Tjn+l _ jjn / pr _ Pvr ^ 

U i,j,k - U i,j,k ^i+l/2,j,ft *i-l/2,j,k) 

\ r i.] + l/2,k r i,j-l/2,k) 



■ A/Q,.,.,. (17) 

The index n represents the time level, and the time and 
space discretization intervals are indicated by At and 
Ari, AO, and Atp for the r-, 9-, and tp-direction, respec- 
tively. The numerical fluxes along the three coordinate 
directions, F r , F e , and F v . are computed by means of 
Marquina's flux formula 63] . A family of local Riemann 
problems is set up at every cell-interface, whose jumps are 
minimized with the use of a piecewise parabolic recon- 
struction procedure (PPM) which provides third-order 
accuracy in space. 

We note that Godunov-type schemes have also been 
implemented recently in 2D and 3D Cartesian codes de- 
signed to solve the coupled system of the Einstein and 
hydrodynamic equations, as reported in |49l |5C1 l53ll64j . 

D. Elliptic solvers 

In the following we present the three different ap- 
proaches we have implemented in our code to numerically 
solve the system of metric equations (|llf) . We compare 
the properties of these solvers with special focus on issues 
like 

- radius and order of convergence, 

- scaling with resolution in various coordinate directions, 

- imposition of boundary conditions, 

- assumptions about the radial extension of the grid, 

- computational performance, 

- parallelization issues, and 

- extensibility from two to three spatial dimensions. 

In order to formalize the metric equations we define a 
vector of unknowns 

u = u p = (<j), acj), fi 1 ^ 2 ,^). (18) 

Then the metric equations (|llf) can be written as 

f(u) = f q (u p ) = 0, (19) 

with / = f q denoting the vector of the five metric equa- 
tions for u (p, q = 1, . . . , 5). For metric solvers 1 and 2 the 
metric equations are discretized at cell centers (r^, 9j, Lp k ) 
on the finite difference grid. Correspondingly, for met- 
ric solver 3 the metric equations are evaluated at col- 
location points (fi,9j,tpk) on the spectral grid. Thus, 
when discretized, Eq. i|19|) transforms into the follow- 
ing coupled nonlinear system of equations of dimension 
5 x n r x ng x n v or 5 x n r x hg x h v , respectively: 

f(u) = faki&mn) = fUM,m,n) = 0, (20) 

with the vector of discretized equations / = fi,j,k = ffj u 
for the unknowns u = U;, m ,„ = u f m „■ F° r this system 
we have to find the roots. Note that, in general, each dis- 
cretized metric equation /? ,, couples both to the other 
metric equations through the five unknowns (indices p), 
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and to other (neighboring) cell locations on the grid (in- 
dices I, m, n). 

All three metric solvers are based on iterative methods, 
where the new value for the metric u s+1 is computed 
from the value at the current iteration s by adding an 
increment Ait 5 which is weighted with a relaxation factor 
/ r . The tolerance measure we use to control convergence 
of the iteration is the maximum increment of the solution 
vector on the grid the iteration is executed on, i.e. 

A< ax = max (Au») = max (A<j k ). (21) 

1. Multidimensional Newton-Raphson solver (Solver 1) 

Solver 1, which was already introduced in the core 
collapse simulations reported in |lll Il2f . uses a multi- 
dimensional Newton-Raphson iteration method to find 
the roots of Eq. i|20|) . Thus, solving the nonlinear system 
is reduced to finding the solution of a linear problem of 
the same dimension during each iteration. The matrix 
A defining the linear problem consists of the Jacobi ma- 
trix of / and additional contributions originating from 
boundary and symmetry conditions (see pd| for further 
details). As the spatial derivatives in the metric equa- 
tions (which also contain mixed derivatives of second or- 
der) are approximated by second-order central differences 
with a three-point stencil, A has a band structure with 
1 + 2d 2 bands of blocks of size 5x5, where d is the 
number of spatial dimensions of the finite difference grid. 
Furthermore, matrix A is sparse and usually diagonally 
dominated. 

A simple estimate already shows that the size n x n 
of the linear problem grows impractically large in 3D. A 
resolution of 100 grid points in each coordinate direction 
results in a square (5x 10 6 ) x (5x 10 6 ) matrix A. Thus, di- 
rect (exact) inversion methods, like Gauss-Jordan elimi- 
nation or exact LU decomposition, are beyond practical 
applicability, as these are roughly n 3 processes, where n 
is the dimension of the matrix. Even when exploiting 
the sparsity and band structure of A the linear problem 
remains too large to be solved on present-day comput- 
ers in a reasonable time by using iterative methods like 
successive over-relaxation (SOR) or conjugate gradient 
(CG) methods with appropriate preconditioning. 

Because of these computational restrictions, the use of 
solver 1 is restricted to 2D axisymmetric configurations, 
where the matrix A has nine bands of blocks. Even in 
this case, for coupled spacetime and hydrodynamic evolu- 
tions, the choice of linear solver methods is limited: The 
computational time spent by the metric solver should not 
exceed the time needed for one hydrodynamical time step 
by an excessive amount. We have found that a recursive 
block tridiagonal sweeping method [6^ (for the actual 
numerical implementation, see [TT| ^1 yields the best per- 
formance for the linear problem. Here the three leftmost, 
middle, and rightmost bands are combined into three new 
bands of n r blocks of size (5 x ng) x (5 x ng) and which are 



inverted in a forward-backward recursion along the bands 
using a standard LU decomposition scheme for dense ma- 
trices. Actual execution times for this method and the 
scaling with grid resolution are given in Section TlVB II 

We point out that the recursion method provides 
us with a non-iterative linear solver, and the Newton- 
Raphson method exhibits in general very rapid and ro- 
bust convergence. Therefore, solver 1 converges rapidly 
to an accurate solution of the metric equations (|19|l even 
for strongly gravitating, distorted configurations, irre- 
spective of the relative strength of the "hydrodynamics" 
term Sh and "metric" term S m in the metric equations 
(see Eq. <|12f> ). Its convergence radius is sufficiently large, 
so that even the flat Minkowski metric can be used as an 
initial guess for the iteration, and the relaxation factor 
f x can be set equal to 1. Note that in solver 1 every 
metric function is treated numerically in an equal way; 
in particular, the equations for each of the three vector 
components of the shift vector j3 l are solved separately. 

In its current implementation, solver 1 exhibits a par- 
ticular disadvantage, which will be discussed in more de- 
tail in Section TlV B 21 As its spatial grid, on which the 
metric equations are discretized, is not radially compact- 
ified, there is a need for explicit boundary conditions of 
the metric functions u at the outer radial boundary of 
the finite difference grid. This poses a severe problem, as 
there exists no general analytic solution for the vacuum 
spacetime surrounding an arbitrary rotating fluid config- 
uration in any coordinate system. Even in spherical sym- 
metry, our choice of isotropic coordinates yields equations 
with noncompact support terms, which leads to imprecise 
boundary conditions, as demonstrated in Section III B 31 
Therefore, as an approximate boundary condition for an 
arbitrary matter configuration with gravitational mass 
M g , we use the monopole field for a static TOV solution, 

M i - Ms. 

2r i + 

evaluated at rfd. The influence of this approximation on 
the accuracy of the solution for typical compact stars is 
discussed in Section TlV B 21 We emphasize that the use 
of a noncompactified finite radial grid is not an inherent 
restriction of this solver method. However in the case of 
metric solver 1, for practical reasons we have chosen to 
keep the original grid setup as presented in |llj . where 
both the metric and hydrodynamic equations are solved 
on the same finite difference grid. 

Finally, a further drawback of solver 1 is its inefficiency 
regarding scalability on parallel or vector computer archi- 
tectures. The recursive nature of the linear solver part of 
this method prevents efficient distribution of the numeri- 
cal load onto multiple processors or a vector pipeline. In 
combination with the disadvantageous scaling behavior 
of the linear solver with resolution (see also Table ITTTl be- 
low), these practical constraints render any extension of 
solver 1 to 3D beyond feasibility. 
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2. Conventional iterative integral nonlinear Poisson solver 
(Solver 2) 

While solver 1 makes no particular assumption about 
the form of the (elliptic) equations to be solved, solver 2 
exploits the fact that the metric equations l|ll|) can 
be written in the form of a system of nonlinear cou- 
pled equations with a Laplace operator on the left hand 
side (|12[1 . A common method to solve such kind of equa- 
tions is to keep the right hand side S(u) fixed, solve 
each of the resulting decoupled linear Poisson equations, 
Aii s+1 — S(u s ), and iterate until the convergence crite- 
rion 12]) is fulfilled. 

The linear Poisson equations are transformed into inte- 
gral form by using a three-dimensional Green's function, 

-^//^//sin.w/v«^#, (23) 
At: J J J \x-x'\ 

where the spatial derivatives in S are approximated by 
central finite differences. The volume integral on the 
right hand side of Eq. 12: il) is numerically evaluated by ex- 
panding the denominator into a series of radial functions 
fi(r,r') and associated Legendre polynomials P ; m (cos0), 
which we cut at I — 10. The integration in Eq. 1)23(1 . 
which has to be performed at every grid point, yields a 
problem of numerical size (n r x ng x n v ) 2 . However, the 
problem size can be reduced to n r xng x n v by recursion. 
Thus, solver 2 scales linearly with the grid resolution 
in all spatial dimensions (see Section TlVB II) . However, 
while the numerical solution of an integral equation like 
Eq. i|23|) is well parallelizable, the recursive method which 
we employ to improve the resolution scaling performance 
poses a severe obstacle. In practice only the paralleliza- 
tion across the expansion series index I (or possibly cyclic 
reduction) can be used to distribute the computational 
workload over several processors. 

An advantage of solver 2 is that it does not require 
the imposition of explicit boundary conditions at a finite 
radius due to the integral form of the equations. De- 
manding asymptotic flatness at spatial infinity fixes the 
integration constants in Eq. 1|23[1 . However, as the metric 
equations contain in general source terms with noncom- 
pact support (see Section Hi B 3fl . the radial integration 
must be performed up to infinity to account for the source 
term contributions. As the discretization scheme used in 
solver 2 limits the radial integration to some finite ra- 
dius ny, the metric equations are solved only approxi- 
mately if the source terms with noncompact support are 
nonzero. The consequences of this fact are discussed in 
Section IIV B 21 As in the case of metric solver 1, the 
metric solver 2 could be used with a compactified radial 
coordinate as well. 

One major disadvantage of solver 2 is its slow con- 
vergence rate and a small convergence radius. For sim- 
plicity, we decompose the metric vector equation for the 
shift vector (3 l into three scalar equations for its compo- 
nents. If the 0-component of the shift vector does not 



vanish, f3 2 ^ 0, and if the spacetime is nonaxisymmetric, 
solver 2 does not converge at all (probably due to diverg- 
ing terms like (3 / sin 2 9 in the vector Laplace operator). 
Even when using a known solution obtained with another 
metric solver as initial guess, solver 2 fails to converge. 
Thus, the use of solver 2 is limited to axisymmetry. Even 
so, when /3 2 ^ 0, a quite small relaxation factor / r ss 0.05 
is required. Furthermore, as the iteration scheme is of 
fix-point type, it already has a much lower convergence 
rate than e.g. a Newton-Raphson scheme. Both factors 
result in typically several hundred iterations until conver- 
gence is reached (see Section llVB 1|) . For strong gravity, 
the small convergence radius restricts the initial guess to 
a metric close to the actual solution of the discretized 
equations. 



3. Iterative spectral nonlinear Poisson solver (Solver 3) 



The basic principles of this iterative solver are similar 
to the ones used for solver 2: A numerical solution of the 
nonlinear elliptic system of the metric differential equa- 
tions is obtained by solving the associated linear Pois- 
son equations with a fix-point iteration procedure until 
convergence. However, instead of using finite difference 
scalar Poisson solvers, solver 3 is built from routines of 
the publicly available Lorene library an d uses spec- 
tral methods to solve scalar and vector Poisson equa- 
tions |6q . 

Before every computation of the spacetime metric, the 
hydrodynamic and metric fields are interpolated from the 
finite difference to the spectral grid by the methods de- 
tailed in Section UlI B 21 All three-dimensional functions 
are decomposed into Chebyshev polynomials T n (r) and 
spherical harmonics Y l m (9, ip) in each domain. When us- 
ing solver 3 the metric equations JSJ) are rewritten in 
order to gain accuracy according to the following trans- 
formations. The scalar metric functions 4> and a have 
the same type of asymptotic behavior near spatial infin- 
ity, ^Ir-.oo — 1 + A^(r), a| r _oo ~ 1 + Aa(r), with A^(r) 
and Aa(r) approaching as r — > oo. Therefore, to ob- 
tain a more precise numerical description of the (usually 
small) deviations of <f> and a from unity, we solve the 
equations for the logarithm of (f> and a<p, imposing that 
lnc/i and ln(a</>) approach zero at spatial infinity. Another 
important difference to the other two solvers is that the 
vector Poisson equation for the shift vector (j % is not de- 
composed into single scalar components, but instead the 
entire linear vector Poisson equation is solved, including 
the jV'Vfc operator on the left hand side. Therefore, the 
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system of metric equation to be solved reads 

/ K 
A\n<j ) =-4ircl ) 4 [phW 2 -P> 10 



V" 16tt 
-V 1 \n<j) Viln^, 



/ IK- K^' 

A In a0 = 2^0 4 p/i(3VF 2 - 2) + 5P + --^ ) (24) 

\ lo7r 

-V 1 lna0 Vi lnc«/>, 
A/3* + \v l V k (3 k = lGna^S* + 2,r"A"V / " 



J ^6 



During each iteration a spectral representation of the 
solution of the linear scalar and vector Poisson equa- 
tions associated with the above system is obtained. The 
Laplace operator is inverted (i.e. the linear Poisson equa- 
tion is solved) in the following way: For a given pair of 
indices I and m of Y l m (9, ip), the linear scalar Poisson 
equation reduces to an ordinary differential equation in 
r. The action of the differential operator 



dr 2 



2d_ 

r dr 



1(1 + 1) 



(25) 



acting thus on each multipolar component (/ and m) of 
a scalar function corresponds to a matrix multiplication 
in the Chebyshev coefficient space. The corresponding 
matrix is inverted to obtain a particular solution in each 
domain, which is then combined with homogeneous solu- 
tions (r and 1/r , for a given I) to satisfy regularity and 
boundary conditions. The matrix has a small size (about 
30 x 30) and can be put into a banded form, owing to 
the properties of the Chebyshev polynomials, which facil- 
itates its fast inversion. For more details about this pro- 
cedure, and how the vector Poisson equation is treated, 
the interested reader is addressed to 66]. Note also that 
when solving the shift vector equation, /3 1 is decomposed 
into Cartesian components defined on the spherical polar 
grid (see |ofij). 

The spatial differentials in the source terms on the 
right hand sides of the metric equations are approximated 
by second-order central differences in solvers 1 and 2, 
while they are obtained by spectral methods in solver 3 
(see Section lTlIB 1|) . When using ~ 30 collocation points, 
very high precision (~ 10~ 13 ) can be achieved in the eval- 
uation of these derivatives. Another advantage of metric 
solver 3 is that a compactified radial coordinate u = 1/r 
enables us to solve for the entire space, and to impose ex- 
act boundary conditions at spatial infinity, u = 0. This 
ensures both asymptotic flatness and fully accounts for 
the effects of the source terms in the metric equations 
with noncompact support. Solver 3 uses the same fix- 
point iteration method as solver 2, but does not suf- 
fer from the convergence problem encountered with that 
solver. Due to the direct solution of the vector Poisson 
equation for the shift vector (3 Z , it converges to the cor- 
rect solution in all investigated models (including highly 
distorted 3D matter configurations with velocity pertur- 
bations, see Section flVB 1(1 . Furthermore, this can be 



achieved with the maximum possible relaxation factor, 
/ r = 1, starting from the flat metric as initial guess. 

However, the strongest reason in favor of solver 3 is 
its straightforward extension to 3D. As mentioned previ- 
ously, both metric solvers 1 and 2 are limited to axisym- 
metric situations. The spectral elliptic solvers provided 
by the Lorene library are already intrinsically three- 
dimensional. Indeed, even in axisymmetry the spectral 
grid of solver 3 requires h v = 4 grid points in the in- 
direction order to correctly represent the Cartesian com- 
ponents of the shift vector. 

There is an additional computational overhead due to 
the communication between the finite difference and the 
spectral grids. These computational costs may actually 
become a dominant part when calculating the metric 
(as will be shown in Section IIV A(l . The interpolation 
methods also have to be chosen carefully to obtain the 
desired accuracy. Furthermore, spectral methods may 
suffer from Gibbs phenomena if the source terms of the 
Poisson-like equations contain discontinuities. For the 
particular type of simulations we are aiming at, disconti- 
nuities are present (supernova shock front, discontinuity 
at the transition from the stellar matter distribution to 
the artificial atmosphere at the boundary of the star). 
This can result in high-frequency spurious oscillations of 
the metric solution, if too few radial domains are used, or 
if the boundaries of the spectral domains are not chosen 
properly. As mentioned before, a simple way to reduce 
the oscillations is to filter out part of the high-frequency 
spectral coefficients. 

As the C++ routines of the Lorene library in the 
current release are optimized for neither vector nor par- 
allel computers, solver 3 cannot yet exploit these archi- 
tectures. However, we were able to improve the computa- 
tional performance by coarse-grain parallelizing the rou- 
tines which interpolate the metric solution in the spectral 
representation to the finite difference grid. 



E. Extraction of gravitational waves 

In a conformally fiat spacetime the dynamical grav- 
itational wave degrees of freedom are not present |ll| . 
Therefore, in order to extract information regarding the 
gravitational radiation emitted in core collapse events 
and in rotating neutron star evolutions, we have imple- 
mented in the code the 3D generalization of the axisym- 
metric Newtonian quadrupole formula used in |l0l ITU 

El- 

Note that we use spherical polar components for the 
tensors of the radiation field. 

Whereas in axisymmetry there exists only one inde- 
pendent component of the quadrupole gravitational ra- 
diation field hjj- in the transverse traceless gauge, 



hl T (r 7 9) = -A + (0)e 



+ • 



(26) 
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in three dimensions we have 



hJ7 



(r,6,<p) = ^[A + (9,cp)e++A x (9,<p)e x ], (27) 



with the unit vectors e+ and e x dehned as 



e> 



ee 



+ e u 



ee- 



(28) 
(29) 



The amplitudes A + and A x are linear combinations 
of the second time derivative of some components of the 
quadrupole moment tensor /y, which for simplicity we 
evaluate at tp = on the polar axis and in the equatorial 
plane, respectively: 



^4+ — hi — I22, 
A p x =2J 12 , 

A c + = I33 — I22 , 
A c x =-27i 3 , 



at 9 = (pole), 



(30) 



at 9 = 7r/2 (equator). (31) 



A direct numerical calculation of the quadrupole mo- 
ment in the standard quadrupole formulation, 



ho = / dVp 



is, 



(32) 



results in high frequency noise completely dominating the 
wave signal due to the presence of the second time deriva- 
tives in Eq. Ij31|l . Therefore, we make use of the time- 
differentiated quadrupole moment in the first moment of 
momentum density formulation, 



L = / dv P * 



-Sij (vixi 



+ v 2 x 2 + V3X3) 



and stress formulation, 



(33) 

dVp* [2viVj - ad,® - Xjdi®], (34) 



of the quadrupole formula |67l l68| . 

In the above equations, x% and Vi are the coordinates 
and velocities in Cartesian coordinates, respectively. 
When evaluating Eq. (|34l) numerically, we transform Vi to 
spherical polar coordinates. In the quadrupole moment, 
we use p* = pW(j) e instead of p as in [l(J, HH H2] , as this 
quantity is evolved by the continuity equation (note that 
both quantities have the same Newtonian limit). This 
also allows a direct comparison with the results presented 
in which we show in Section llVB4l For a discussion 
about the ambiguities arising from the spatial derivatives 
of the Newtonian potential <& in Eq. (|34|l in a general 
relativistic framework and their solution (which we also 
employ in this work), we refer to |T^ . 

The total energy emitted by gravitational waves can 
be expressed either as a time integral, 



2 

E gw = — I at 



-hxhz — hih3 — I22I33 



+h\ 



/si 



:>Ui2 



hi 



I23 



(35) 



or, equivalently, as a frequency integral, 



1 
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v 2 dv 



'I11I22 — I11I33 — I22I; 



'3:? 



-2 ~2 .-.2 /.-.2 ~2 .-.2 

-In + I22 + I33 + 3 J12 + ^13 + ^23 



, (36) 



where hj{v) is the Fourier transform of hj{t). We point 
out that the above general expressions reduce to the fol- 
lowing ones in axisymmetry: 

A\ = 0, A p x = 0, A% = I, Al = 0, (37) 
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E gw = — I dtl = — I v*dvl 



15 



(38) 



with I = I33 — I22 being the only nonzero independent 

.-.2 

component of the quadrupole tensor, and / being the 
• ■2 

Fourier transform of / . The quadrupole wave amplitude 
Af 2 used in [1113, El is related to I according to A^ = 

We have tested the equivalence between the wave- 
forms obtained by the axisymmetric code presented 
in 0,0,^3 and those by the current three-dimensional 
code using the corresponding axisymmetric model. In all 
investigated cases, they agree with excellent precision. 



IV. CODE TESTS AND APPLICATIONS 

We turn now to an assessment of the numerical code 
with a variety of tests and applications. We recall that we 
do not attempt in the present paper to investigate any re- 
alistic astrophysical scenario, which is deferred to subse- 
quent publications. Instead, we focus here on discussing 
standard tests for general relativistic three-dimensional 
hydrodynamics code, which were all passed by our code. 
In particular, we show that the code exhibits long-term 
stability when evolving strongly gravitating systems like 
rotational core collapse and equilibrium configurations of 
(highly perturbed) rotating relativistic stars. Each sepa- 
rate constituent methods of the code (HRSC schemes for 
the hydrodynamics equations and elliptic solvers based 
on spectral methods for the gravitational field equations) 
has already been thoroughly tested and successfully ap- 
plied in the past (see e.g. [44L l45l l66| and references 
therein). Therefore, we mainly demonstrate here that 
the coupled numerical schemes work together as desired. 



A. Interpolation efficiency and accuracy 

The interpolation procedure from the finite difference 
grid to the spectral grid has been described in Sec- 
tion IIII B 21 Among the three possible algorithms we 
have implemented in the code, the most efficient turned 
out to be the one based on a piecewise parabolic in- 
terpolation (see Table . It is as fast as the piecewise 



14 



TABLE I: Execution time ifd^ sp and accuracy A/fd^ sp for 
the interpolation of a test function / t (r, 8, ip) (see text) from 
the finite difference grid to the spectral grid, listed for different 
finite difference grid resolutions n r xngxn^ and interpolation 
types. The interpolation methods are piecewise linear (type 
1), piecewise parabolic (type 2), and globally minimizing the 
norm of the second derivative of the interpolated function |4§1| 
(type 0). The spectral grid has a resolution of h r = 17, hg — 
17, and h v — 16 grid points. 



n r X ng X n v 


Type 


_i_ r 1 
ifd^sp [sj 


A/fd-» sp 


Lq norm] 


400 x 200 x 800 


2 


5.13 


5.0 


X 


10" 8 


400 x 200 x 800 


1 


5.12 


7.0 


X 


10" 6 


400 x 200 x 800 





9.44 


1.8 


x 


1Q ~6 


400 x 200 x 400 


2 


2.92 


3.1 


x 


lO" 7 


400 x 200 x 200 


2 


1.43 


1.6 


x 


KT 6 


400 x 200 x 100 


2 


0.77 


1.7 


x 


10" 5 


400 x 200 x 10 


2 


0.09 


1.3 


x 


10" 2 


400 x 100 x 800 


2 


2.55 


3.1 


x 


lO" 7 


400 x 50 x 800 


2 


1.60 


1.8 


X 


10" 6 


400 x 5 x 800 


2 


0.32 


2.0 


X 


10~ 3 


200 x 200 x 800 


2 


3.61 


2.7 


X 


lO" 7 


100 x 200 x 800 


2 


1.81 


2.1 


X 


10" 6 


50 x 200 x 800 


2 


1.40 


1.6 


X 


10" 5 


5 x 200 x 800 


2 


0.99 


1.4 


X 


10~ 2 



linear interpolation, and more accurate than the algo- 
rithm based on the minimization of the second deriva- 
tive of the interpolated function. Table [i] shows, for 
a particular example of an interpolated test function 
/t(r, 9, <p) = exp [— r 2 (l 4- sin 2 (9 cos 2 ip)] , the relative ac- 
curacy A/i n t (in the L norm) achieved by this interpo- 
lation, as well as the CPU time spent on a Pentium IV 
Xeon processor at 2.2 GHz. The spectral grid consists of 
two domains (nucleus + shell) with h r = 17, hg = 17, 
and h v = 16. The outer radius of the nucleus is located 
at 0.5, and the outer boundary of the shell is at 1.5 (cor- 
responding to the radius of the finite difference grid rfd). 

This test demonstrates that the piecewise parabolic 
interpolation is indeed third-order accurate, and that the 
time spent scales roughly linearly with the number of 
points of the finite difference grid in any direction. We 
have made other tests which show that the interpolation 
accuracy is independent of n, and that it scales in time 
like O (n 3 ) + O (n 3 ) , where h and n are the number of 
points used in each dimension by the spectral and the 
finite difference grid, respectively. The interpolation is 
exact, up to machine precision, for functions which can 
be expressed as polynomials of degree < 2 with respect 
to all three coordinates. 

The direct spectral summation from the spectral to 
the finite difference grid is a very precise way of evalu- 
ating a function: For smooth functions, the relative er- 
ror decreases like exp(— h) (infinite order scheme). This 
property is fulfilled in our code, as shown in Table ITT1 for 
the same test function /t(r, 9, if) and the same domain 



TABLE II: Execution time t sp ^fd and accuracy A/ sp ^fd for 
the evaluation of a test function f t (r,8,<j>) (see text) on the 
finite difference grid from its representation in spectral co- 
efficients, listed for different numbers of spectral grid points 
h r x he x h v . The finite difference grid has a resolution of 
n r — 100, ng — 50, and n v — 30 grid points. 



h r x hg x h v 


isp^fd [s] 


A/sp^fd [Lq norm] 


33 x 17 x 64 


75.8 


1.5 x 10~ 15 


33 x 17 x 32 


38.4 


5.5 x 10~ 9 


33 x 17 x 16 


19.6 


2.6 x 10" 4 


33 x 17 x 8 


10.3 


2.8 x 10~ 2 


33 x 9 x 64 


40.8 


6.4 x 10~ 9 


33 x 5 x 64 


23.4 


3.2 x 10" 4 


17 x 17 x 64 


41.2 


1.9 x 10 -13 


9 x 17 x 64 


24.6 


9.2 x 10~ 7 


5 x 17 x 64 


16.7 


1.9 x 10" 3 



setup as for Table |U (again the timings are for a Pen- 
tium IV Xeon processor at 2.2 GHz). Double precision 
accuracy is reached with a reasonable number of points 
(h r = 33, hg — 17, and h v = 64). According to Table UTI 
the CPU cost scales linearly with the number of coeffi- 
cients h in any direction. We have also confirmed that 
it scales linearly with the number of finite difference grid 
points n in any direction. The drawback of this most 
straightforward procedure is that it requires O (h 3 n 3 ) 
operations, which is much more expensive than the in- 
terpolation from the finite difference grid to the spectral 
one, and even more expensive than the iterative proce- 
dure providing the solution of system (|2"4"|> . Nevertheless, 
it is computationally not prohibitive since the overall ac- 
curacy of the code does not depend on n (which can thus 
remain small) . A way to reduce the execution time is to 
use a partial summation algorithm (see e.g. |62j|). which 
needs only O (nn 3 ) + O (n 2 n 2 ) + O (n 3 n) operations, at 
the additional cost of increased central memory require- 
ment. Another alternative is to truncate the spectral 
sum, staying at an accuracy level comparable to that of 
finite difference differential operators. 



B. Solver comparison in 2D 

1. Convergence properties 

The theoretical considerations about the convergence 
properties of the three implemented metric solvers (as 
outlined in Section IIII Dp are checked by solving the 
spacetime metric for a 2D axisymmetric rotating neu- 
tron star model in equilibrium (labeled model RNS), 
which we have constructed with the method described 
in Komatsu et al. [7(j • This model has a central density 
p c = 7.905 x 10 14 g cm" 3 , obeys a polytropic EoS with 
7 = 2 and K — 1.455 x 10 5 (in cgs units), and rotates 
rigidly at the mass shedding limit, which corresponds to 
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FIG. 2: Comparison of the convergence behavior for the three 
metric solvers in 2D. For solver 1 (filled circles), the maximum 
increment A«J,„ per iteration s decreases to the threshold 
Au" hr = 10~ 15 (lower horizontal dotted line) within less than 
10 iterations, while solver 3 (asterisks) needs more than 40 
iterations to reach its (less restrictive) threshold (upper hor- 
izontal dotted line) of 10 -6 . The very low relaxation factor 
needed for solver 2 (filled squares) results in a remarkably slow 
convergence, requiring more than 700 iterations. The solid 
lines mark the approximate linear decrease of log Auf nax . 



a polar-to-equatorial axis ratio of 0.65. These model pa- 
rameters are equivalent to those used for neutron star 
models in H3,|7ll 

To the initial equilibrium model we add an r- and 9- 
dependent density and velocity perturbation, 



p = Pmi 1 + 0.02 sin 2 (n^j (l + sin 2 (26<)) 
v r = 0.05 sin 2 hr^ (l + sin 2 (26»)) , 
v g = 0.05 sin 2 (n^j sin 2 (2<9), 
v v = V v ini + 0.05 sin 2 (rr^j (l + sin 2 (20)) , 



(39) 



where r s is the ((9-dependent) stellar radius, and v r — 
\J Viv 1 , vg = \fv2V 2 , and v v = V v^v 3 . The metric equa- 
tions (Eqs. ltTT|) for solvers 1 and 2, and Eqs. (|2~4"|) in the 
case of solver 3) are then solved using the three imple- 
mented metric solvers. The perturbation of v r and vg 
ensures that the metric equations yield the general case 
of a shift vector with three nonzero components, which 
cannot be obtained with an initial model in equilibrium. 

We point out that by adding the perturbations speci- 
fied in Eq. 1)39|1 and calculating the metric for these per- 
turbed initial data, we add a small inconsistency to the 
initial value problem. As the Lorentz factor W in the 
right hand sides of the metric equations contains met- 
ric contributions (which are needed for computing the 
covariant velocity components), it would have to be it- 



erated with the metric solution until convergence. How- 
ever, as the perturbation amplitude is small, and as we 
do not evolve the perturbed initial data, we neglect this 
small inconsistency. 

The most relevant quantity related to convergence 
properties of the metric solver is the maximum increment 
Au^ ax of all metric components on the grid (see Fig- - 
As expected solver 1 exhibits the typical quadratic de- 
cline of a Newton-Raphson solver to its threshold value 
Au* hr = 10~ 15 . As the methods implemented in solvers 2 
and 3 correspond to a fix-point iteration, the decline of 
their metric increment is significantly slower. Therefore, 
for the Poisson-based solvers, we typically use a less re- 
strictive threshold Auj* hr = 10~ 6 . While the spectral 
Poisson solver 3 allows for a relaxation factor of 1 and 
thus for a still quite rapid convergence, the conventional 
Poisson solver 2 requires more than 700 iterations due 
to its much smaller relaxation factor imposed by the (3 2 - 
equation. 

It is worth stressing that all three solvers show rather- 
robust convergence, if one keeps in mind that the initial 
guess is the flat spacetime metric. If the metric is chang- 
ing dynamically during an evolution, the metric values 
from the previous computation can be used as new start- 
ing values, which reduces the number of iterations by 
about a factor of two with respect to those reported in 
Fig. El 

Besides the convergence rate, the execution time t m 
required for a single metric computation and its depen- 
dence on the grid resolution is also of paramount rel- 
evance for the practical usefulness of a solver. These 
times for one metric computation of the perturbed RNS 
stellar model on a finite difference grid with various r- 
and (9-resolutions on an IBM RS/6000 Power4 proces- 
sor are summarized in Table IIHI As theoretically ex- 
pected, both solver 1 and 2 show a linear scaling of t m 
with the number of radial grid points n r , i.e. the ratio 
fn r = t m (n r ) /t m (n r / 2 ) is approximately 2. While the 
integration method of solver 2 shows linear dependence 
also for the number of meridional grid zones ng, the in- 
version of the dense ng x ng matrices during the radial 
sweeps in solver 1 is roughly a ng process. Thus, the the- 
oretical value of r ng — 8 for that solver is well met by the 
results shown in Table IITTI We note that for even larger 
values of ng, specific processor properties like cache- miss 
problems can even worsen the already cubic scaling of 
solver 1, while for ng > 64 solver 2 fails to converge 
altogether. On the other hand for solver 3 t m is approx- 
imately independent of the number of finite difference 
grid points in either coordinate direction, as the number 
of spectral collocation points is fixed. A dependence on 
n r and ng can only enter via the interpolation procedure 
between the two grids, the time for which is, however, 
entirely negligible in 2D. 

The break even point for the three solvers corresponds 
roughly to a resolution of 100 x 32 grid points at t m ~ 
20 s. We emphasize that this value of t m is much larger 
than the time needed for one hydrodynamic step at the 



TABLE III: Metric solver execution time t m for different fi- 
nite difference grid resolutions n r x ng for the three metric 
solvers in 2D applied to the perturbed rotating neutron star 
model RNS. The ratios a nr (a ng ) between execution times for 
a given n r (ng) and for half that resolution exhibit the be- 
havior expected from theoretical considerations. The spectral 
grid has a resolution of h r — 33, ng = 17, and h v = 4 grid 
points. 





Solver 1 




Solver 2 




Solver 3 




n r x ng 


tm [s] 


a nr 


o, ng 


tm. [s] 


fln r 




t m [s] 


a„ r 




50 x 16 


1.8 






2.8 






20.7 






100 x 16 


3.7 


2.0 




5.9 


2.1 




20.6 


1.0 




200 x 16 


7.4 


2.0 




12.9 


2.2 




20.8 


1.0 




50 x 32 


12.5 




6.9 


5.9 




2.1 


20.8 




1.0 


100 x 32 


25.4 


2.0 


6.9 


12.3 


2.1 


2.1 


20.5 


1.0 


1.0 


200 x 32 


50.8 


2.0 


6.9 


27.1 


2.2 


2.1 


21.7 


1.1 


1.0 


50 x 64 


109.7 




8.8 


12.4 




2.1 


20.9 




1.0 


100 x 64 


224.2 


2.0 


8.8 








21.5 


1.0 


1.1 


200 x 64 


445.2 


2.0 


8.8 








21.7 


1.0 


1.1 



same resolution, which is roughly th ~ 0.1 s. From the 
results reported in Table UTTl it becomes evident that due 
to the independence of t m on the finite difference grid 
resolution in the spectral metric solver 3, this method 
is far superior to the other two solvers for simulations 
requiring a large number of grid points in general, and 
particularly in (9-direction. 

2. Radial fall-off of the metric components 

When comparing in Section AlI Dl the theoretical foun- 
dations of the three alternative metric solvers imple- 
mented in the code, we already raised the issue of the 
existence of source terms with noncompact support in 
the metric equations I|11J) (see Section lll B 3|) . Neither the 
Newton-Raphson-based solver 1, which requires explicit 
boundary conditions at the finite radius rfd (which are in 
general not exactly known and possibly time-dependent) , 
nor the conventional iterative Poisson solver 2, which in- 
tegrates the Poisson-like metric equations only up to the 
same finite radius rfd, are able to fully account for the 
nonlinear source terms, even if the radial boundary of the 
finite difference grid is in the vacuum region outside the 
star, rfd > r s . 

Hence, both solvers yield a numerical solution of the 
exact metric equations only in very few trivial cases, like 
e.g. the solution for the metric of a spherically symmet- 
ric static matter distribution (TOV solution), when the 
metric equations reduce to Poisson-like equations with 
compact support. However, due to the radial compacti- 
fication of the spectral grid, which allows for the Poisson 
equations to be numerically integrated out to spacelike 
infinity, the spectral solver 3 can consistently handle all 
noncompact support source terms in the metric equa- 
tions in a non-approximative way. This property holds 
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FIG. 3: Equatorial profile of the shift vector component /3 vo 
obtained by different metric solvers compared with the correct 
profile from the initial data solver (solid line) for the rotating 
neutron star model RNS. Due to its approximate boundary 
value, the profile from solver 1 (dashed line) shows large devi- 
ation from the correct solution, particularly for a grid bound- 
ary rfd close to the stellar equatorial radius r s e (upper panel) . 
As solver 2 (dashed-dotted line) needs no explicit boundary 
conditions, its solution matches well with the correct solution, 
with improving agreement as rfd is at larger distance from r s e 
(lower panel) . The compactified radial grid of solver 3 (dotted 
line) fully accounts for non-compact support terms, and thus 
agrees very well with the correct solution, independent of the 
location of rfd. The radii r sc and rfd are indicated by vertical 
dotted lines. 



even when the metric quantities are mapped from the 
spectral grid onto the finite difference grid, the latter ex- 
tending only to rfd. Thus, we expect that only solver 3 
captures the correct radial fall-off behavior of the metric 
quantities outside the matter distribution. 

In the following we illustrate the effects of noncompact 
support terms in the metric equations on the numerical 
solution using the three different solvers. Fig. [3] shows 
the radial equatorial profiles of the rotational shift vector 
component (3 V = ^/733/3 3 for the rapidly rotating neutron 
star initial model (RNS) specified in Section TlVB II ob- 
tained with the three alternative metric solvers. While 
we restrict our discussion to the particular metric quan- 
tity /3tp e we notice that the radial fall-off behavior and 
the dependence on the solver method is equivalent for all 
other metric components. 

In the upper panel of Fig. |3| the equatorial stellar 
boundary r sc is very close to the radial outer boundary 
of the finite difference grid, r se = 0.9 rfd (both indicated 
by vertical dotted lines). The star and the exterior at- 
mosphere are resolved using n rs = 90 radial grid points 
for the star and n ra = 10 radial grid points for the at- 
mosphere (along the equator), respectively, and rig = 30 
meridional points. The spectral solver 3 uses n r = 33 



17 



radial and hg — 17 meridional grid points. 

If the boundary value for the metric at rfd is exact, 
solver 1 always yields the correct solution, irrespective of 
the source terms not having compact support. For sta- 
tionary solutions like rotating neutron stars these exact 
values can in principle be provided by the initial data 
solver. However, for instance in a dynamical situation, 
exact values cannot be provided, and we are forced to use 
approximate boundary conditions, which we choose ac- 
cording to Eq. (|22J) . As the approximate boundary value 
for solver 1, (3 v (r{d) = 0, is far from the exact value, 
the corresponding profile of the shift vector (dashed line) 
strongly deviates from the correct /3 ve obtained by the 
initial data solver (solid line). Note that the exact solu- 
tion is given only for r < r so , due to limitations of the 
initial solver method [7(j • As shown in the lower panel of 
the figure, with increasing distance of the finite difference 
grid boundary from the stellar boundary (rfd = 2.0 r sc 
with n rs = n ra — 90), the approximation for /3 ve (7"fd) 
improves noticeably, and so does the matching of (3 ve 
with the correct solution. 

On the other hand, as the integral approach of solvers 2 
and 3 requires no specific boundary conditions at a fi- 
nite radius (contrary to solver 1), the numerical solution 
for (3 ve agrees well with the correct solution even for an 
integration boundary rfd close to the stellar boundary 
r sc (dashed-dotted and dotted lines in Fig. |3 respec- 
tively). For rfd S> r sc , when the influence of the source 
terms with noncompact support is increasingly picked up 
by the radial integral, the solutions supplied by solver 2 
rapidly approach the correct one. The terms with non- 
compact support usually do not contribute strongly to 
the solution of the metric equations (except in cases of 
very strong gravity and extremely rapid contraction or 
rotation). Thus, solver 2 is superior to solver 1 when 
approximate boundary values must be used, Eq. (|22|l . 
Solver 3, on the other hand, has the key advantage over 
solver 2 of using very accurate spectral methods for solv- 
ing the Poisson equation over the entire spatial volume 
due to its compactificd radial coordinate. Hence, irre- 
spective of the distance of rfd from r so , it yields the same 
results on the finite difference grid, onto which the results 
are mapped from the spectral grid. 

The (small) difference between the results for (3 V c from 
solver 3 and from the initial data solver is partly due to 
the accuracy of the numerical schemes and the mapping 
between different grids, and particularly a result of the 
CFC approximation of the field equations employed by 
the evolution code (note that the initial data are gen- 
erated from a numerical solution of the exact Einstein 
metric equations). In the case of rapidly rotating neu- 
tron star models we have found that the truncation error 
and the error arising from the mapping of the initial data 
to the evolution code is typically more than one order of 
magnitude smaller than the error which can be attributed 
to the CFC approximation, if a grid with a resolution 
n r ~ 100, ng ~ 30 and n r = 33, fig — 17 is used. For es- 
timates of the quality of the CFC approximation in such 



cases, see [llj and references therein. 

We again note that, in principle, the use of a compact- 
ified radial grid is not confined to the spectral solver 3. 
A finite difference grid extending to spatial infinity could 
be used for solvers 1 and 2 as well. However, in that 
case either the exterior atmosphere would also have to 
be extended to the entire grid too (generating unneces- 
sary computations), or only the relevant portion of the 
grid containing the star would have to be evolved in time 
(creating an additional boundary). When using solver 3, 
there is a clearcut split between the finite difference grid 
and the spectral grid. Thus, the hydrodynamic quanti- 
ties can be defined on a grid with an atmosphere of only 
small size, while the metric in the compactified domain 
can be computed very accurately with only few radial 
collocation points due to the exponential convergence of 
spectral methods in this smooth region. Additionally, 
the Lorene library provides the use of a compactified 
radial domain as an already implemented option at no 
extra cost. 



3. Axisymmetric core collapse to a neutron star - 
Construction of the spectral grid domains 

As all three metric solvers yield equally precise nu- 
merical solutions of the spacetime metric in 2D, they 
give nearly identical results when applied to simula- 
tions of rotational core collapse, as shown in Fig. 0] 
For the results presented in this figure we have chosen 
the stellar core collapse model labeled A3B2G4 in [l^ 
(model SCC in the following), which rotates differen- 
tially and moderately fast, and has an initial central 
density p c — 10 10 g cm -3 . The initial adiabatic index 
is reduced from y t = 4/3 to 71 = 1.3 during contrac- 
tion, and is increased to 72 = 2.5 beyond supranuclear 
matter densities, p > p uuc — 2.0 x 10 14 g cm -3 . The 
details of the EoS for this model are given by Eq. JBJ. 
As the metric calculation is computationally very ex- 
pensive, it is done only every 100/10/50 hydrodynamic 
time steps before/during/after core bounce, and extrapo- 
lated in between (for details on the satisfactory accuracy 
of this procedure see ^j). The number of zones used 
in the finite difference grid is n r — 200 and ng = 30, 
with logarithmic spacing in the r-direction and a central 
resolution of 500 m, and an equidistant spacing in the 
^-direction. Again, the grid resolution of the spectral 
solver 3 is n r = 33 and ng = 17. 

In the upper panel of Fig. 0] we plot the time evolution 
of the central conformal factor C , which rises steeply 
when the central density increases to supranuclear den- 
sities, reaches a maximum at the time of core bounce t\, 
(vertical dotted line) , and subsequently approaches a new 
equilibrium value with decreasing ringdown oscillations. 
This new state, which is reached asymptotically, signals 
the formation of a pulsating compact remnant which can 
be identified with the nascent proto-neutron star. Each of 
the three curves in this upper panel is the result of using 
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FIG. 4: Time evolution of the central conformal factor <f> c 
(upper panel) for the core collapse model SCC, using metric 
solver 1 (solid line), 2 (dashed line), and 3 (dashed-dotted 
line), respectively. All three solvers yield similar results. The 
small relative differences of less than 10~ 3 in cf> c (lower panel) 
obtained with solvers 1 and 3 (solid line) and solver 2 and 3 
(dashed line) prove that numerical variations of the metric 
from each solver are of the order of the small overall dis- 
cretization error of the entire evolution code. The time of 
bounce is indicated by the vertical dotted line. 



one of the three available metric solver (see caption for 
details). The lower panel of the figure demonstrates that 
the relative differences found in the dynamical evolution 
of our representative core collapse model are negligibly 
small when using either metric solver, which proves the 
applicability of any of the metric solvers in 2D. 

However, in such a highly dynamical situation, where 
the relevant radial scales vary by a factor of about 100, 
solver 3 requires a special treatment of the radial do- 
main setup of the spectral grid defined in Section lill B II 
During the infall phase of a core collapse simulation the 
contracting core must be suffiently resolved by the ra- 
dial grid, and thus we adjust the radius of the nucleus 
rd dynamically before core bounce. (Note that this is no 
contradiction to the assumption /<j = const, in Eq. (|15fl . 
as /d may change between subsequent metric calculations 
during the evolution.) Initially the value of rd is given 
by half the stellar radius. As the evolution proceeds it is 
set equal to the radial location of the sonic point in the 
equatorial plane (once unambiguously detected). Alter- 
natively r<i can be determined by the radius enclosing a 
shell of a fixed fraction of the total rest mass of the star 
(typically 10%), whereby rd moves inward during the col- 
lapse, too. In either case rd is held fixed when some mini- 
mal radial threshold rd m in is crossed, which we set equal 
to the radius of some given radial grid point (e.g. the 
40th grid point at r^). This ensures that there is always 
a sufficient number of grid points on the finite difference 
grid, such that the interpolation to the spectral grid is 
well behaved. For rid — 6 domains, both approaches 




t [ms] 



FIG. 5: Two different methods for determining the domain 
radii of the spectral grid boundary. The upper panel shows 
the time evolution of the domain radius parameter rd for the 
core collapse model SCC, where rd is either set by the sonic 
point method (solid line; sonic point first detected at t ~ 
23 ms) or by the rest mass fraction method (dashed line). The 
boundary of the finite difference grid rfd , the stellar equatorial 
radius r sc , the minimal domain radius rdmin (set to r4o), and 
the approximate location of shock formation r s h are indicated 
by horizontal dotted lines. The relative difference between the 
values of <f> c from simulations using the two methods (lower 
panel) is less than 10 -4 throughout the evolution. The time 
of bounce tb is indicated by the vertical dotted line. 



yield equally accurate results, the relative difference be- 
tween the values of <fi c being less than 10 -4 throughout 
the evolution of the collapse model SCC (see lower panel 
of Fig. EJ. 

At least for core collapse simulations, the appropriate 
choice of the radial spectral domain setup parameters 
and r<i(t) is crucial, as exemplified in Fig. The re- 
duction of rd with time must follow the contraction of 
the core to a sufficiently small radius, while rdmin must 
retain enough grid points for the nucleus. Furthermore, 
when splitting the spectral grid into several radial do- 
mains, well-behaved differential operators (in particular, 
the Poisson operator) are only obtained if, for a shell- 
type domain, the criterion of thin shell- type domains, 
/d < 2, is fulfilled. This restriction for the ratio /d be- 
tween the outer and the inner radii originates from the 
requirement to keep the condition number of the matrix 
representing (for a given multipolar momentum I) the ra- 
dial Poisson operator (|25p. which is a very fast growing 
function of /d, lower than ~ 10 3 . 

In particular Fig. shows that if rd is not properly 
adjusted or if rd m i n is too large, the central conformal 
factor deviates strongly from the correct value (upper 
panel) . In addition, if the number of domains is too small 
while keeping the radial resolution n r = 33 fixed, the 
conformal factor inside the core shows large amplitude 
oscillations after core bounce, due to a too large value 
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FIG. 6: Importance of the correct spectral domain setup for 
highly dynamic simulations, shown for the core collapse model 
SCC. If the domain radius parameter is not reasonably ad- 
justed (upper panel), e.g. r^ is held fixed at 10% of the initial 
stellar equatorial radius (dashed line), or if the minimal do- 
main radius is too large, ramta — r ioo (dashed-dotted line), 
the central conformal factor <f> c deviates strongly from the cor- 
rect value (solid line; c.f. Fig. 2J. If the number of domains 
is too small (lower panel), e.g. rid = 3 (dashed line) instead 
of rid — 6 (solid line), the metric inside the star (here the 
equatorial conformal factor (f>r wa e at the 100th radial grid 
point) shows strong oscillations after core bounce. The time 
of bounce tb is indicated by the vertical dotted line. 



of /d (lower panel). If /d < 2 is violated because of too 
few domains in a collapse situation, such oscillations are 
even present if the radial resolution n r is increased. 

On the other hand, in quasi-stationary situations with 
no large dynamical radial range (e.g. oscillations of neu- 
tron stars), one can safely reduce rid from 6 to 3 and keep 
rd fixed throughout the evolution. The optimal number 
of domains is thus determined by balancing radial res- 
olution and the requirement of thin shell-type domains 
against computational costs. 



4- Axisymmetric core collapse to a neutron star - 
Comparison with fully general relativistic simulations 

Only recently, fully general relativistic simulations of 
axisymmetric rotational core collapse have become avail- 
able 0] . We now estimate the quality of the CFC ap- 
proximation adopted in our code by s imulating one of the 
core collapse models presented in 17] and comparing the 
results. 

In their simulations, Shibata and Sekiguchi |17| make 
use of the Cartoon method jT^I which reduces the di- 
mensionality of a code based on 3D Cartesian coordinates 
to 2D in the case of axisymmetric configurations. Using 
this approach, and solving the full set of BSSN metric 
equations, these authors present a series of rotational core 
collapse models with parameters close (but not exactly 



equal) to the ones simulated by Dimmelmeier et al. |12fl . 
As an additional difference, p* = pW<j) 6 is employed 
by [17j in the gravitational wave extraction with the first 
moment of momentum density formula, while in the 
wave extraction is performed with the stress formula us- 
ing the density p (see Section Hll El for details). Further- 
more, in the simulations reported in [TtI ]. the equidistant 
Cartesian finite difference grid is repeatedly remapped 
during the collapse, so that the grid spacing in the cen- 
ter increases from initially ~ 3 km to ~ 300 m during 
core bounce. As the outer boundary moves in accord- 
ingly, matter leaves the computational grid, resulting in 
a mass loss of about 3%. 

In their paper, Shibata and Sekiguchi investigated a 
core collapse model which is identical to our model SCC 
(A3B2G4 in [12j) with the exception of a slightly smaller 
rotation length parameter A = A/r sc = 0.25 (compared 
to A = 0.32 in 0]) in the initial equilibrium model. 
They found that the evolution of this model (labeled 
SCCss hereafter) computed with their fully general rela- 
tivisitc code agrees qualitatively well with the evolution 
of our model SCC simulated with our CFC code. How- 
ever, it produces an increased gravitational wave ampli- 
tude of about 20% at the peak during core bounce, and 
up to a factor 2 in the ringdown. Furthermore, the damp- 
ing time of the ringdown signal of model SCCss as shown 
in 17] is significantly longer compared to that of model 
SCC presented in [l|. 

Shibata and Sekiguchi offer several possible explana- 
tions for this noticeable disagreement, the most plausi- 
ble ones being the different functional forms of the rest 
mass density used in the wave extraction method, and 
the different formulations (stress formulation H34|) versus 
first moment of momentum density formulation l|33|) ). 
By comparing waveforms obtained from evolutions of os- 
cillating neutron stars (as presented in |69j), both using 
the quadrupole formula and by directly reading off met- 
ric components, they find that the quadrupole formula 
underestimates the wave amplitude of model SCCss by 
~ 10%. Extrapolating these results they arrive at the es- 
timate that the waveforms presented in |l2j are accurate 
at best to within ~ 20%. Shibata and Sekiguchi claim 
that other differences, namely the CFC approximation 
versus the BSSN formulation, different grid setups, co- 
ordinate choices and slicing conditions, or the small dis- 
crepancy of A in the initial model, have only negligible 
impact on the waveform. 

To test this conjecture, we have simulated the evo- 
lution of model SCC with our new version of the CFC 
code in 2D, and extracted the wave amplitude A c + using 
the first moment of momentum density formulation l|33(l 
with p, and also alternatively substituting p by p* . As 
our results show (see upper panel of Fig. 0), the use of 
p* results in a small increase of A e + by about 20% dur- 
ing the bounce and the ringdown phase, limiting pos- 
sible deviations due to the difference in the quadrupole 
formula stated in ^3 to about 20%. However, the re- 
sults depicted in Fig. [7\ exclude that the doubling of A e + 
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FIG. 8: Influence of differences in the initial model on the 
evolution of the central density p c for rotational core collapse. 
Changing from the collapse model SCC (solid line) to SCCss 
(dashed line) only slightly shifts the time of bounce tb (indi- 
cated by the vertical dotted line), but leads to much stronger 
post-bounce ring down oscillations. Nuclear matter density 
p nuc is indicated by the horizontal dotted line. 



FIG. 7: Influence of the density used in the wave extrac- 
tion equations (upper panel) and of small differences in the 
initial model (lower panel) on the gravitational waveforms 
from rotational core collapse. If p* = phW 2 is used in the 
quadrupole formula (solid line) instead of p (dashed line), the 
wave amplitude A c + increases by about 20% at core bounce 
(upper panel). A change from model SCC (solid line) to model 
SCCss (dashed line), which corresponds solely to a difference 
in the initial configuration, results in a qualitatively differ- 
ent waveform, in particular during the ringdown phase (lower 
panel) . The times of bounce tb are indicated by the vertical 
dotted lines. 

observed by 01 for the ringdown signal is due to the 
wave extraction method. On the contrary, comparing the 
waveforms for model SCC and SCCss (see lower panel of 
Fig. UJ, both computed with our CFC method, shows 
that the strong qualitative difference found by Shibata 
and Sekiguchi is clearly due to the differences in the core 
collapse initial model, notably the small decrease of the 
differential rotation length scale A in model SCCss- This 
gives rise to an approximately 50% higher peak value of 
the amplitude during bounce, and a strong increase of the 
post bounce wave amplitude, as also observed by Shibata 
and Sekiguchi (compare with Fig. 13 (b) in |17|). 

Furthermore, from the evolution of the central density 
computed with our code (see Fig. 0, it is evident that 
model SCCss exhibits significantly stronger ringdown os- 
cillations than model SCC with a somewhat longer damp- 
ing timescale, which is also in good agreement with the 
results in [Tg (see their Fig. 7 (b)). Clearly the small dif- 
ference in the rotation length parameter A of the initial 
model has a major impact on the post-bounce dynam- 
ics of the dense core, which is in turn reflected in the 
gravitational wave signal. 

We have also simulated the evolution of models SCC 
and SCCss using a larger number of radial and merid- 
ional grid points (n r = 250 and ng — 60 with a cen- 
tral radial resolution A rc = 250 m) as compared to 
the standard grid setup with n r = 200, ng = 30, and 



A rc = 500 m (in either case the spectral grid resolution is 
h r = 33 and hg = 17). Neither improving the resolution 
of the finite difference grid nor discarding a significant 
mass fraction in the outer parts of the star (to mimic the 
mass loss introduced by the regridding method in ^tJ) 
have a significant impact on the collapse dynamics or the 
waveform for both initial models. When simulating the 
same collapse model, the observed small differences to 
Shibata and Sekiguchi's results in e.g. the central den- 
sity or the waveform are most likely due to the use of the 
CFC approximation for the spacetime metric employed 
in our code. Nevertheless, for core collapse simulations, 
the results obtained using either CFC or the full Einstein 
equations agree remarkably well. 

C. Applications of the spectral solver 3 in 3D 

1. Computation of a nonaxisymmetric spacetime metric 

While the previous tests were all restricted to 2D (and 
thus solvers 1 and 2 could as well be used), the genuine 
3D properties of the spectral metric solver 3 can be fully 
exploited and tested when applied to the computation of 
the metric for a nonaxisymmetric configuration. For this 
purpose we consider now the uniformly rotating neutron 
star initial model RNS (see Section TlVB 1|) to which we 
add a nonaxisymmetric perturbation. This is done by 
generalizing the expressions in Eq. i|39|) through the mul- 
tiplication of a (/^-dependent term of the form (1 + sin 2 ip). 
The effect of such a perturbation on representative quan- 
tities is depicted in Fig. EI The metric equations (|24|l are 
then integrated using solver 3. Convergence is reached af- 
ter about 50 iterations (threshold value Auj hr = 10 -6 ), 
and the solution for the metric is interpolated from the 
spectral to the finite difference grid. 

To exclude convergence to an incorrect solution and 
errors within the interpolation routine, we compare the 
left and right hand sides, lhs„ and rhs u , of selected metric 
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FIG. 9: Nonaxisymmetric density and velocity perturbation 
of the rapidly rotating neutron star equilibrium model RNS. 
By applying the perturbations described in the text, the origi- 
nal profiles (dashed lines) of the density p along the azimuthal 
direction ip (upper panel), the radial velocity v r along the 
meridional direction 6 (center panel), and the rotation ve- 
locity v v along the radial direction r (lower panel) become 
strongly distorted (solid lines). The 99-dependence of p in 
the upper panel shows the nonaxisymmetric character of the 
perturbation. 



components u on the finite difference grid, in Fig. ^| We 
note that in this figure, along each of the profile direc- 
tions, the two other coordinates are kept fixed [r = ?'5o, 
6 = 7r/4, and <p = 0, respectively). The left and right 
hand sides of the metric equations 124|) for the conformal 
factor tfi and the shift vector components (3 1 and /3 3 , when 
evaluated on the finite difference grid, match very accu- 
rately along all three coordinate directions. The largest 
deviations are found near the rotation axis (0 = 0) for 
P. 

The accuracy of the metric calculation can be better 
quantified by plotting the relative difference of the left 
and right hand sides, A ro i u = |lhs u /rhs u — 1|, rather than 
lhs„ and rhs u alone. This is shown for the metric quan- 
tities <fi, P 1 , and /3 3 in the insets of Fig. Along any of 
the plotted profiles, the spectral solver yields a solution 
for which the relative difference measure is better than 
10~ 2 . As lhs u and rhs u contain second spatial derivatives 
of the metric, evaluated by finite differencing, this is an 
accurate numerical result. We note that some of the met- 
ric components are close to zero or change sign. Hence, 
the relative difference may become large or develop a pole 
at some locations, as can be seen in the insets of Fig.llUI 

Under idealized conditions (i.e. without discontinuities 





FIG. 10: Left (solid line) and right (dashed line) hand sides 
(computed on the finite difference grid) of the equation for the 
metric components cf> along the azimuthal direction ip (upper 
panel), (3 1 along the meridional direction 6 (center panel), 
and /3 3 along the radial direction (lower panel). Even for 
strong nonaxisymmetric perturbations of the rotating neutron 
star model RNS, the metric solver 3 yields a highly accurate 
matching, such that the lines almost lie on top of one another. 
The insets show the relative difference A rc i u between the left 
and right hand sides of the equation for the same metric com- 
ponents. The relative differences are < 1CF 2 , except where 
they exhibit a pole. 



in the source terms of the metric equations, no artifi- 
cial atmosphere, only laminar matter flows, uniform grid 
spacing of the finite difference grid, and perturbations 
which are regular at the grid boundaries), such a test 
case also offers an opportunity to examine the order of 
convergence of the metric solver 3 on the spectral and fi- 
nite difference grid, respectively. To this end we perform 
a metric calculation using increasingly finer resolutions 
on the two grids. By varying the number of spectral col- 
location points in all three spatial directions while keep- 
ing the number of finite difference grid points fixed (at 
high resolution), we observe an exponential decrease of 
the relative differences A re i u between the left and right 
hand sides of the equation for the various metric com- 
ponents u. Correspondingly, the metric solution evalu- 
ated on the finite difference grid exhibits second order 
convergence with grid resolution for a fixed (and high) 
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spectral grid resolution. Furthermore, the (at least) sec- 
ond order accurate time integration scheme of the code 
in combination with the PPM reconstruction of the Rie- 
mann solver also guarantees second order convergence 
during time evolution. For fixed time steps we actually 
observe this theoretical convergence order globally and 
even locally (except close to the grid boundaries, where 
symmetry conditions and ghost zone extrapolation spoil 
local convergence). 

In the three-dimensional case the computational load 
of the interpolation from the spectral grid to the finite 
difference grid after every metric calculation on the spec- 
tral grid becomes significant. The time spent in the in- 
terpolation between grids can, in fact, even surpass the 
computational costs of the spectral metric solution itself 
(see Section HV A|) . As a consequence, the independence 
of the metric execution time t m on the number of finite 
difference grid points found in the axisymmetric case (as 
shown in Table ITTTf) cannot be maintained. Table HVI re- 
ports the summary of runtime results for a single met- 
ric computation of the above neutron star model on an 
IBM RS/6000 Power4 processor. These results indicate 
an (albeit sublincar) increase of t m with the number of 
finite difference grid points. As expected, a doubling of 
the spectral grid resolution e.g. in the (/^-direction (while 
keeping fi r = 33 and ng — 17 fixed) results in a pro- 
portional increase of t m . The runtime scaling results re- 
ported in Table IIVI also demonstrate that the different 
coordinate directions contribute equally to the computa- 
tional costs. 

It is worth pointing out that the other two metric 
solvers we have available in the code fail to compute the 
metric for the nonaxisymmetric neutron star configura- 
tion considered in this section due to the known limita- 
tions (excessive computing time for solver 1, convergence 
problems for solver 2). 



Stability of symmetric configurations against 
perturbations 



TABLE IV: Dependence of the metric solver execution time 
t m on the finite difference grid resolution n r xng xn v and the 
spectral grid azimuthal resolution h v using the metric solver 3 
in 3D for the nonaxisymmetrically perturbed rotating neutron 
star model RNS. For typical finite difference grid point num- 
bers, the ratio r ntp between execution times for a given n v 
and for half that resolution is smaller than 2, i.e. the increase 
of tm is less than linear. Furthermore, when doubling both 
the radial and meridional grid zones, a sublinear increase in 
the corresponding ratio r„ r e < 4 is observed. Doubling the 
spectral resolution increases t m by ~ 2. For com- 
parison, the values of tm for the corresponding axisymmetric 
model are given at the bottom. 
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model RNS. Model SNS has the same central density and 
EoS as model RNS described in Section llVB II To each 
equilibrium model SNS and RNS we respectively add an 
axisymmnetric (r, 6)- and a nonaxisymmetric (r, in- 
dependent three-velocity perturbation of the form 



V r = 0.02 Sh^ 7T 



(1 + a sin 2 (2(9)) 



vg = 0.02 sin 2 (tt— ) asm 2 (260, 



(40) 



An important requirement for any hydrodynamics 
code is the preservation of the symmetry of an initially 
symmetric configuration during time evolution. In a 
practical application this means that if a small pertur- 
bation is added to symmetric and stable initial data, the 
perturbation amplitude must not grow in time. Due to 
the choice of spherical polar coordinates (r,0,ip), our 
code is particularly well suited to test the preservation 
of the symmetry of spherically symmetric and axisym- 
metric initial data. Additionally, this coordinate choice 
implies that when simulating axisymmetric or spherically 
symmetric problems, either one or two dimensions can be 
trivially suppressed, respectively, which results in consid- 
erable savings of computational time. 

Next, we present results from the evolution of both 
a spherically symmetric neutron star model (labeled 
SNS) and the axisymmetric rapidly rotating neutron star 



and 



v r = 0.02 sin 



(n^j (1 + sin 2 (26>)) (l + a sin 2 ip) 



(41) 



0.02 sin 2 tt— (l + sin 2 (26»))(l + asin 2 <p) 



v g = 0.02 sin 2 ( tt— ) sin 2 (26») (l + a sin 2 ip) , 



V<p — V <p ini ~\~ ' 



respectively, where a is the perturbation amplitude. 
Model SNS is then evolved in time using the code in 
axisymmetric 2D mode, and model RNS using the fully 
3D capabilities of the code. The metric is calculated ev- 
ery 100 (300) time steps in 2D (3D) and extrapolated in 
between. The number of finite difference grid zones is 
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12 in the 3D case. Correspondingly, for 
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FIG. 11: Time evolution of a symmetry violating perturba- 
tion. The upper two panels correspond to the spherically 
symmetric model SNS, and the lower two panels to the ax- 
isymmetric model RNS. The relative variation in density Ap 
(solid line), radial velocity Av r (dashed line), rotational ve- 
locity Av v (dotted line), and conformal factor Acj> (dashed- 
dotted line) show a remarkable constancy in time (note that 
Av v is nonzero only for the rotating model RNS). The sym- 
metry violating variation of the different fields scale with the 
initial perturbation amplitude (horizontal dotted lines; left 
panels: a — 10 -3 ; right panels: a = 10 -6 ). 



the spectral grid we use n r = 25, ng — 13, n v = 4 in 2D, 
and n r = 25, ng = 13, h v = 6 in 3D. 

The results of the evolution of the symmetry violating 
perturbations in both models are depicted in Fig. 1111 The 
upper panels correspond to model SNS which is evolved 
up to 5 ms, while the bottom panels correspond to model 
RNS which is only evolved up to 1 ms. The left and right 
panels differ by the value of the initial amplitude a of 
the velocity perturbation. We observe that the pertur- 
bation amplitude, measured as the relative difference Aq 
of an arbitrary matter or metric quantity q evaluated at 
two points of constant r (for model SNS) and constant 
r, 8 (for model RNS), remains practically unchanged for 
many hydrodynamic time scales. Note that the spikes 
in Aq appearing in Fig. ^2 are the poles associated with 
a vanishing q. Fig. ^] a ls° shows that the amplitude of 
the symmetry violation Aq approximately scales with the 
amplitude a of the initial velocity perturbation (indicated 
by horizontal dotted lines). 

In the course of many hydrodynamic time scales, the 
perturbations (which are of small amplitude, a <C 1) will 
be finally damped due to the intrinsic numerical viscos- 
ity of the schemes implemented in the code. However, if 
the rotation rate (3 of a rotating neutron star were high 
enough such that /? > p s or even (3 > /3d, perturbations 




FIG. 12: Radial profile of the equatorial rotation velocity v ve 
for the unperturbed axisymmetric rapidly rotating neutron 
star model RNS evolved in 3D. The profile of v vc at t = 10 ms 
(dashed line) closely reproduces the initial profile (solid line). 
The stellar equatorial radius r sc and the boundary of the finite 
difference grid r^ are indicated by vertical dotted lines. 



of the form given by Eq. Ij41|l could trigger the onset of 
physically growing modes, leading to bar mode instabili- 
ties. 



3. Evolution of an axisymmetric uniformly rotating 
neutron star in 3D 

The ability to handle long-term evolutions of rapidly 
rotating relativistic equilibrium configurations is a diffi- 
cult test for any numerical code. To demonstrate the ca- 
pabilities of our code to pass this stringent test we evolve 
the rotating neutron star initial model RNS in 3D until 
t = 10 ms, which corresponds to about 10 hydrodynamic 
time scales and rotation periods. The simulation is per- 
formed with a resolution for the finite difference grid of 
n r = 100, ng = 30, n v = 8, and h r — 33, fig = 17, n v = 6 
for the spectral grid. During the evolution, the metric is 
calculated every 100 time steps and extrapolated in be- 
tween. 

The preservation of the radial profile of the rotation ve- 
locity Vtp e along the equator over a long evolution time is 
shown in Fig. 1121 Depicted is the initial equilibrium solu- 
tion (solid line) as a function of the radial coordinate (in 
the equatorial plane) and the final configuration (dashed 
line), after an evolution time of 10 ms (about 10 rota- 
tional periods). The figure shows that v v remains close 
to its initial value in the interior of the star, showing the 
strongest (but still small) deviations near the stellar sur- 
face (at the interface to the artificial atmosphere). This 
local decrease of v v due to interaction of stellar matter 
with the atmosphere and its depencence on the order of 
the reconstruction scheme has also been observed in other 
studies (see e.g. [7l|). 

It is important to emphasize that the accurate preser- 
vation of the rotational profile is achieved because of 
the use of third-order cell-reconstruction schemes for the 
hydrodynamics equations, such as PPM, as first shown 
by UM ■ Despite the comparably coarse resolution of the 
finite difference grid and the use of the CFC approxima- 
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FIG. 13: Time evolution of the radial velocity at half the stel- 
lar equatorial radius v r he for the perturbed rapidly rotating 
neutron star model RNS. The radial velocity shows regular 
oscillations with neither a noticeable drift nor damping when 
the 3D code is used in low resolution (solid line) as well as 
for the 2D code with high resolution (dashed line). For com- 
parison, the dashed-dotted line shows v r he when no explicit 
perturbation is added. In this case the oscillations are trig- 
gered by truncation errors and (mostly) by the error resulting 
from using the CFC approximation in the evolution code. 

tion for the gravitational field equations, our code cap- 
tures the profile of v vc at the stellar boundary about 
as accurately as codes solving the full Einstein metric 
equations coupled to the hydrodynamics equations |50| . 
or codes restricted to hydrodynamic evolutions in a fixed 
curved spacetime (i.e. using the so-called Cowling ap- 
proximation) |7l| . 

Long-term evolutions of rotating neutron stars as the 
one presented here can be effectively used for extract- 
ing the oscillation frequencies of the various pulsation 
eigenmodes of the star. This topic has been traditionally 
studied using perturbation theory (see e.g. |73[ and refer- 
ences therein). In recent years fully nonlinear hydrody- 
namical codes have helped to drive progress in the field. 
They have provided the quasi-radial mode-frequencies of 
rapidly rotating relativistic stars, both uniformly and dif- 
ferentially rotating, which is a problem still not amenable 
to perturbation techniques (see e.g. |50L iTll ITl. I75ll7r| ) . 

In order to test our code against existing results we 
show next an example of the procedure to compute mode- 
frequencies using the model RNS. The frequencies can in 
principle be extracted from a Fourier transform of the 
time evolution of various pulsating quantities when the 
oscillations are triggered by numerical truncation errors. 
However, the results significantly improve when a pertur- 
bation of some specific parity is added to the initial equi- 
librium model. To excite small amplitude quasi-radial 
oscillations, we hence apply an I = radial velocity per- 
turbation to the equilibrium configuration of the form 

v r = asm 2 (^—^j i ( 42 ) 

with an amplitude a = —0.01. 

Due to this perturbation, various metric and hydro- 
dynamic quantities exhibit very regular periodic oscilla- 
tions around their equilibrium state, as shown for the 



TABLE V: Comparison of the oscillation frequencies of two 
perturbed equilibrium neutron star models SNS and RNS 
with different axis ratios r sp /r sc obtained with the current 
code (both in 2D and 3D) and with the Cactus code |5(3| . 
The frequencies for the fundamental mode /p and for the first 
harmonic /hi computed with the current code show a relative 
difference with respect to the Cactus code (in parentheses) 
of at most 2%. Due to the coarse spatial resolution used, the 
3D code results were only calculated to 3 significant figures. 





SNS 


RNS 




^sp/^sc 


= 1.00 




= 0.65 


Code 


/f [kHz] 


fm [kHz] 


/f [kHz] 


/hi [kHz] 


current (3D) 


1.40 (3.4) 


3.95 (0.2) 


1.20 (0.4) 


3.68 (1.0) 


current (2D) 


1.463(0.9) 


3.951(0.2) 


1.219(2.0) 


3.659 (1.6) 


Cactus 


1.450 


3.958 


1.195 


3.717 



radial velocity v r in Fig. 1131 The pulsations, which show 
no noticeable numerical damping during the entire dura- 
tion of the simulation (10 ms), are extracted at half the 
stellar equatorial radius. The same oscillation pattern is 
obtained when instead of using the 3D code (solid line in 
the figure) the model is evolved using the code in axisym- 
metric mode (dashed line in Fig. ll3l with finite difference 
grid size of n r = 160, ng — 60). The latter, axisym- 
metric setup is currently being used in a comprehensive 
parameter study of the oscillation frequencies of rotating 
neutron star models Note that Fig. ^| also demon- 
strates that the oscillation amplitude scales linearly with 
the initial perturbation amplitude a (at least if a <C 1), 
which was chosen as a = —0.005 in the 2D simulations. 
In the radial velocity, neither an offset nor a noticeable 
drift with time can be observed. This is in agreement 
with previous results using alternative formulations and 
different numerical codes |50l l7l| . 

Time evolution data like the one shown in Fig. 1131 
can be used to extract the eigenmode frequencies. A 
Fourier transformation of different metric and hydrody- 
namic quantities at various locations in the star yields 
identical (discrete) frequencies. Table IVl summarizes the 
frequencies f-p and /hi for the quasi-radial fundamental 
mode and its first harmonic overtone, respectively. Both 
frequencies obtained with the current 3D code differ only 
by a few percent from those computed with the code in 
2D or the Cactus code, which is based on a Carte- 
sian grid and uses the BSSN formulation for the Einstein 
equations [5(i| . 

Additionally, we have investigated the influence of grid 
resolution and finite evolution time on the accuracy of the 
frequency extraction. We have found that the differences 
in the frequencies between the 2D and 3D simulations 
presented in Table E] can be almost entirely attributed 
to the twice as long evolution time of the 2D simulation 
(20 ms), for which the Fourier transformation renders 
more accurate frequencies. For practical evolution times 
of several tens of milliseconds and for grid resolutions 
better than n r ~ 100 and ng ~ 30, the extracted oscil- 
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lation frequencies are almost independent of the number 
of grid points used. 

Note also that the mode-frequencies agree well even 
though we have used different perturbation amplitudes 
a in the 3D and 2D simulations (while in the Cactus 
run an I = rest mass density perturbation with an 
amplitude a = 0.02 was used). Table Ivl hence proves 
that our code is able to simulate rotating neutron stars 
in a fully three-dimensional context for sufficiently long 
time scales to successfully extract oscillation frequencies. 



4- Evolution of a nonaxisymmetric uniformly rotating 
neutron star in 3D 



Contrary to the small amplitude nonaxisymmetric per- 
turbations employed in Section IIV C 21 we turn now to 
assess the ability of the numerical code to manage long- 
term stable evolutions of strongly gravitating systems 
with large departures from axisymmetry. This is an es- 
sential test for future astrophysical applications of the 
code as e.g. the numerical investigation of bar mode in- 
stabilities in rotating neutron stars. 

For this purpose we construct a uniformly rotating 
neutron star model with the same parameters as model 
RNS, but with only half the central density. The fi- 
nite difference grid extends out to rfd = 80 km, with 60 
equidistant radial grid points resolving the neutron star 
out to r sc = 18.6 km. The atmosphere is covered by 80 
logarithmically spaced radial grid points. The number of 
angular zones used in the finite difference grid is ng = 24 
and n v — 32, respectively, while the spectral grid has 
h r = 17, fig = 13, and h v — 12 grid points in 3 radial 
domains. 

On top of the equilibrium neutron star model we add 
a strongly nonaxisymmetric (i.e. (^-dependent) perturba- 
tion of the rest-mass density 



P = Pini + a p c sin 



r 

2r~, 



sin 10 ip for r < 2r s 



(43) 

with an amplitude a = 0.1, which yields an / = m = 2 
bar-like structure. The rotation velocity of the uniformly 
rotating unperturbed neutron star is extrapolated into 
the areas filled with matter by the perturbation. The ini- 
tial configuration with the perturbation added is shown 
in the left panel of Fig. ^] 

We have chosen this particular (albeit unphysically 
strong) perturbation and velocity field in order to pre- 
vent both, an immediate accretion of the added matter 
bars on to the neutron star or an ejection. This allows us 
to follow the rotation of the neutron star for a time com- 
parable to its rotation period (which is about 1 ms for the 
unperturbed neutron star). The density and rotation ve- 
locity plots in Fig.^Jafter t = 0.5 ms (center panel) and 
£ = 1.0 ms (right panel) prove this property of the cho- 
sen perturbation. These plots also demonstrate that the 



corotating bar structures slowly disappear. The inner- 
most parts are being gradually accreted by the neutron 
star, which leads to a significant initial rise in the central 
density, as shown in Fig. 1151 At later times the more mas- 
sive neutron star oscillates with a period of i osc ~ 1.0 ms 
around a new quasi-equilibrium state, which possesses a 
central density of more than 50% above the initial equi- 
librium central density. Despite this strong interaction of 
the bar perturbation with the neutron star, the rotation 
profile inside the neutron star remains uniform through- 
out the evolution, although the rotation velocity nearly 
doubles during the oscillation maxima. This behavior 
is most likely due to the particular choice of a uniform 
rotation profile for the initial bar perturbation. 

For the outer parts of the initial bar, the increasing 
distance from the neutron star and the sufficiently high 
specific angular momentum prevents their accretion onto 
the neutron star. Thus the matter in this region of the 
bar drifts to larger radii during the evolution. As on the 
dynamical timescales considered of one rotation period 
there is no efficient transport mechanism of local angular 
momentum by viscous effects (which act on much longer 
timescales), the evolution leads to the development of 
spiral arms which are clearly visible in the middle and 
right panels of Fig. The outer parts of these arms 
are centrifugally expelled from the finite difference grid, 
crossing the outer boundary at t ~ 0.84 ms. By the end 
of the simulation, at t — 4 ms, there is neither significant 
backscattering of matter from the outermost boundary 
of the radial grid, nor there are numerical artifacts vis- 
ible at the leading or trailing edges of the spiral arms. 
This proves that our numerical treatment of the radial 
boundary conditions and of the artificial low density at- 
mosphere surrounding the star have the desired behavior. 

Fig. E| shows that already after an evolution time of 
~ 1 ms, the evolution of the spiral arms has no further 
significant impact in the dynamics of the neutron star, as 
then the slowly decaying oscillation around the final equi- 
librium state exhibits a rather regular ring-down pattern. 
Plotted in this figure is also the time evolution of the cen- 
tral density for a model with an amplitude a = 0.01 of the 
initial perturbation given by Eq. I)43|) (dashed line). In 
addition, the dashed-dotted line shows the correspond- 
ing time evolution of p c for an unperturbed model (the 
small amplitude oscillations are in this case triggered 
by the truncation errors of the numerical schemes and 
by the use of the CFC approximation in the evolution 
code). The similarity in the behavior of p c in these cases 
demonstrates that for perturbations with an amplitude 
a < 0.01, the dynamics of the central neutron star is 
virtually unaffected by the initial bar and by the spiral 
arms forming at later times. However, we observe that 
also for small values of a spiral arms develop which are 
stable over many rotation periods. 

Apparently, strong nonaxisymmetric perturbations of 
the form (|43|l give rise to significant gravitational wave 
emission. The waveforms of the nonzero gravitational 
wave amplitudes A c +1 A v +1 and A' (as shown in the upper, 
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FIG. 14: Evolution of a strongly distorted nonaxisymmetric rotating neutron star model. The color coded distribution of 
log p on the equatorial plane shows how the initial perturbation (left panel) is partly accreted by the neutron star, and partly 
stretched into spiral arms (center panel). After about one rotation period of the neutron star, the trailing spiral arms have 
grown considerably in size (right panel). The rotation velocity v v is indicated by white arrows. Note that the atmosphere 



(color coded in black) has a density of much less than 
domain are shown. 



10 g cm , and that only the innermost 60 km of the computational 
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FIG. 15: Time evolution of the central density p c for dis- 
torted nonaxisymmetric rotating neutron star models. If the 
distortion is strong (a = 0.1, solid line), matter accretion 
from the rotating bars results in a steep initial increase of p c , 
which slowly settles down to a new equilibrium state (indi- 
cated by the horizontal dotted line). For a small perturba- 
tion (a = 0.01, dashed line), the evolution of p c follows very 
closely that of an unperturbed model (dashed-dotted line). 



center, and lower panel of Fig. 1161 respectively) exhibit 
peak values of up to ~ 15 x 10 3 cm for the model with a 
perturbation amplitude a = 0.1 (solid lines). In Fig. 1161 
we also present the waveforms for the model with a bar 
perturbation of amplitude a — 0.01 (dashed lines). Their 
amplitudes are roughly a factor 10 smaller than those 
of the corresponding waveforms of the model with a = 
0.1. Thus we can infer that the gravitational radiation 
amplitude approximately scales with a. 

We emphasize that owing to the particular form of 
the perturbation (|43|l . the x-mode of the gravitational 



radiation is zero at the equator, A c x — 0. We also note 
that if instead of the nonaxisymmetric perturbation in 
Eq. I|43|) we use an axisymmetric one, 



ap c 



r 

2^ s 



for r < 2r= 



(44) 



then the x-mode of gravitational radiation vanishes com- 
pletely, and only the +- mode is present (dashed-dotted 
line in the upper panel of Fig. 116(1 . Additionally, in 
axisymmetry the +-mode on the pole is always zero, 
0. 

We point out that the waveform pattern for the model 
with the a = 0.1 bar perturbation in Fig. 1161 does not 
solely reflect the oscillation and ring-down structure of 
the central neutron star, as visible in the time evolution 
of p c in Fig. El F° r instance the +- mode at the equator 
(upper panel) decays on a much longer time scale than 
the corresponding ring-down time of p c . On the other 
hand, the waveforms for the two polarizations of the ra- 
diation at the pole exhibit their peaks during the first 
oscillation of p c and then decay rapidly (center and lower 
panel) . However, after an evolution time of ~ 2 ms their 
amplitudes increase again. From this behavior we deduce 
that initially the waveform signal is dominated by the 
gravitational wave emission from the oscillating neutron 
star. As this contribution decays during the ring-down, 
the wave emission from spiral arms becomes increasingly 
important. As they expand into the atmosphere the ra- 
dial weight arm in the quadrupole formula compensates 
for the relatively low density of the spiral arms, and the 
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FIG. 16: Gravitational wave signal for distorted nonaxisym- 
metric rotating neutron star model. If the distortion is strong 
(a — 0.1, solid lines), the nonzero gravitational wave ampli- 
tudes A e + (upper panel), A v + (center panel), and A p x (lower 
panel) reach peak values of up to ~ 15, 000 cm. The ampli- 
tudes reduce significantly for a — 0.01 (dashed lines). If an 
axisymmetric perturbation with a = 0.1 is applied (dashed- 
dotted line), only the ^4+ gravitational wave mode is present. 



radiation emitted in this region becomes visible in the sig- 
nal. We cannot clearly attribute the late-time increase in 
the waveform amplitude to the onset of a bar mode in- 
stability, because the rotation parameter (3 of our model 
clearly falls short of the approximate threshold for dy- 
namical growth of bar modes: (3 ~ 0.14 <C (3^- We plan 
to investigate this issue more thoroughly in the future. 

The maximum amplitude A ~ 15 x 10 3 cm of the wave 
signal for a = 0.1 corresponds to a dimensionless gravi- 
tational wave amplitude /i~5x 10 -19 at a distance of 
r = 10 kpc to the source. Thus, in this case of a strongly 
nonaxisymmetric artificial perturbation, the typical wave 
amplitudes have a value of roughly one order of mag- 
nitude above the ones of waveforms obtained from the 
simplified models of rotational supernova core collapse in 
axisymmctry by Dimmclmcier et al. [12|. For the wave- 
forms plotted in Fig. ^j] we utilize the stress formula l|54"|) 
with p* as density. The use of this formula efficiently 
reduces the numerical noise in the signal as compared 
with the first moment of momentum density formula and 
particularly with the standard quadrupole formula. 

We consider the grid resolution used in this test simu- 
lation to be the minimal one required for obtaining rea- 
sonably converged results. By repeating the same model 
with different grid resolutions we are able to estimate 
that the waveform amplitudes are correctly computed 
within ~ 30% accuracy. 



V. CONCLUSIONS 

In this paper we have presented a new three- 
dimensional general relativistic hydrodynamics code 
which is primarily intended for applications of stellar 
core collapse to a neutron star or a black hole, as well 
as for studies of rapidly rotating relativistic stars which 
may oscillate in their quasi-normal modes of pulsation, 
emitting gravitational radiation, or which may be sub- 
ject to nonaxisymmetric instabilities. The main novelty 
of this code compared to other existing numerical rela- 
tivistic codes is that it combines very accurate state-of- 
the-art numerical methods specifically tailored to solve 
the general relativistic hydrodynamics equations on the 
one hand, and the gravitational field equations on the 
other hand. More precisely, the hydrodynamic equations, 
formulated in conservation form, are solved using high- 
resolution shock-capturing schemes based upon approxi- 
mate Riemann solvers and third-order cell-reconstruction 
interpolation procedures, while the elliptic metric equa- 
tions are solved using an iterative nonlinear solver based 
on spectral methods. Furthermore, the present code also 
departs noticeably from other three-dimensional codes 
in the coordinate system used in the formulation of the 
equations and in the discretization. In our approach both 
the metric and the hydrodynamics equations are formu- 
lated and solved numerically using spherical polar coordi- 
nates. In the present investigation we have adopted the 
so-called conformal flatness approximation of the Ein- 
stein equations, which reduces them to a set of five el- 
liptic nonlinear equations, particularly suited for the use 
of spectral methods. Recently, constrained formulations 
of the full Einstein equations in which elliptic equations 
have a preeminence over hyperbolic equations have been 
reported, and appear to be amenable to the current code. 

The main purpose of the paper has been to assess the 
code by demonstrating that the combination of the finite 
difference grid and the spectral grid, on which the hydro- 
dynamics and metric equations are respectively solved, 
can be successfully accomplished. This approach, which 
we call Mariage des Maillages (French for grid wedding), 
results in high accuracy of the metric solver and, in prac- 
tice, has allowed for fully three-dimensional applications 
using computationally affordable resources, along with 
ensuring long term numerical stability of the evolution. 
To facilitate the Mariage des Maillages, i.e. the combi- 
nation of the finite difference grid for the hydrodynamic 
solver and the spectral grid for the metric solver, a so- 
phisticated interpolation and grid communication scheme 
has been used. In addition, we have compared our novel 
approach to two other, finite difference based, methods 
to solve the metric equations, which we already employed 
in earlier axisymmetric investigations [TTL . 

We have presented a variety of tests in two and three 
dimensions, involving neutron star spacetimes and stellar 
core collapse. Axisymmetric simulations have also been 
performed to compare core collapse to neutron stars us- 
ing the CFC approximation and full general relativity, 
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for which only very recently results have become avail- 
able [T3I . This comparison has shown the suitability of 
the conformally flat approximation for such mildly rela- 
tivistic scenarios. Furthermore, the code has succeeded 
in simulating the highly perturbed nonaxisymmetric con- 
figuration of a uniformly rotating neutron star for several 
dynamical times. This simulation has also been used to 
assess the 3D gravitational waveform extraction capabil- 
ities of the code. In summary the numerical experiments 
reported in the paper demonstrate the ability of the code 
to handle spacetimes with and without symmetries in 
strong gravity. In future work we plan to apply this code 
to simulations of stellar core collapse to neutron stars or 
black holes in three dimensions, and particularly to stud- 
ies of the nonlinear development of bar mode instabilities 
in rapidly rotating neutron stars. 
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APPENDIX A: DIFFERENCES TO PREVIOUS 
2D CFC SIMULATIONS 

1. Compact form of the Euler equation sources 

In the axisymmetric CFC code presented in [TllIT^ the 
source terms Qj for the hydrodynamic momentum equa- 
tions (Euler equations) were evaluated on the finite differ- 
ence grid using a formulation containing time derivatives 
and explicit Christoffel symbols (see Equation J2J): 



Q 3 = T** 1 



dx^ 



r ^9\j 



(Al) 



Using the relation between the Christoffel symbols and 
the derivatives of the spacetime metric, 



a _ 1 xs ( d 9su . dg s ^ dg^ 



\ dx^ dx h 



dx b 



(A2) 



the sources Qj can be written in a more compact form 
as 



1 rp^v ^g^u 

2 dxP ' 



(A3) 



In this formulation, only spatial derivatives of the metric 
are needed, and the numerical evaluation of Qj involves 
significantly fewer terms, making a numerical implemen- 
tation both faster and more accurate. For these reasons, 
we have preferred the use of Equation (|A3|) to Equa- 
tion (|A1|> in the code presented in this paper. 



2. Exact numerical conservation of the 
hydrodynamic equations 

As emphasized in Section 5.4 in 11], the conserved hy- 
drodynamic quantity in the system of conservation equa- 
tions is n °t simply the state vector U but rather y/^yU 
with y/j = cj) 6 r 2 sin 9. Therefore, if only the state vec- 
tor U is evolved, this gives rise to an additional source 
term Q which contains time derivatives of the conformal 
factor <j). These generally time-dependent source terms 
result in a variation of the volume-integrated state vec- 
tor with time, and thus in a violation of exact numerical 
rest mass and angular momentum conservation of sev- 
eral percent, even though the "physical" sources vanish, 
Q = (see Figs. 9 and 10 in Ref. |ll|). 

It is not possible to evolve \FjU in a straightforward 
way and then consistently solve the elliptic metric equa- 
tions 111 1H on the new time slice. This is due to the fact 
that the sources for these equations contain the pressure 
P, which can only be extracted from U but not from 
y/jU. However, one can make use of the time evolution 
equation for the conformal factor, Eq. JJJJ, to obtain an 
auxiliary value for and thus for y/7 on the new time 
slice. With this the state vector U can be consistently 
calculated from y/jU after the time evolution step to the 
new time slice, which in turn is used in the sources of the 
metric equations (| 1 II) . These are subsequently solved on 
the new time slice. With the help of this reformulation of 
the hydrodynamic time evolution problem in the current 
code (in combination with the compact time- independent 
form for the sources in the Euler equations, Eq. (|A3|) ~). 
we are able to achieve exact numerical conservation of 
the total rest mass and angular momentum up to ma- 
chine roundoff errors, provided that there is no artificial 
atmosphere and no mass flow across the outer radial grid 
boundary. 



3. Shift vector boundary conditions 

The results for the evolution of the central density 
p c and the waveform for the core collapse model SCC 
(A3B2G4 in 03) presented in this paper slightly dif- 
fer from those reported in the previous paper by Dim- 
melmeier et al. 12] . This is partly due to the improve- 
ments related to evaluating the Euler equation source 
terms in compact form and using exact numerical conser- 
vation in the new code, as discussed above. However, the 
main reason for the small discrepancy is that in the sim- 
ulations in |12| a symmetric boundary condition for the 
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shift vector component 1 across the equatorial plane was 
chosen. This leads to a nonzero value for (3 2 at 9 = ir/2 
close to and after core bounce, i.e. when meridional mo- 
tions set in. As a consequence of this, the deviation is 
stronger for models where rotation plays a significant role 
in the collapse dynamics. 

The physically accurate antisymmetric equatorial 
boundary condition for 1 which is used in the present 
code, systematically yields lower post-bounce values for 
p c in regular collapse type models compared to the sim- 
ulations presented in jlj, with a difference of 11% on 
average. For models which show multiple bounce behav- 
ior, we obtain a lower p c also at core bounce. 

Accordingly, the waveform amplitudes and frequen- 
cies of the gravitational radiation are altered by a small 



amount (-11% for \Af£ | max and -18% for v). Despite of 
these small quantitative changes, the qualitative state- 
ments related to the influence of general relativistic ef- 
fects in rotational core collapse made by Dimmelmeier et 
al. [12| remain unaffected, even when the antisymmetric 
boundary condition is used. We particularly emphasize 
that the change in the boundary condition for 1 plays 
no role when comparing our results with the fully gen- 
eral relativistic simulations by Shibata and Sekiguchi |l7| 
discussed in Section HV B 41 

We note that for all core collapse models presented in 
the parameter study by Dimmelmeier et al. results 
obtained with the new boundary condition for 1 can be 
found in the revised waveform catalogue pfj . 
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